GitHub Repository: NVIDIA/cuda-q-academic
Path: blob/main/hybrid-workflows/01_quantum_enhanced_optimization_LABS.ipynb
¹¹²⁸ views

Kernel: Python 3 (ipykernel)

Quantum Enhanced Optimization for Radar and Communications Applications

The Low Autocorrelation Binary Sequences (LABS) is an important and challenging optimization problem with applications related to radar, telecommunications, and other signal related applications. This CUDA-Q Academic module will focus on a clever quantum-enhanced hybrid method developed in a collaboration between Kipu Quantum, University of the Basque Country EHU, and NVIDIA for solving the LABS problem.

Other CUDA-Q Academic modules like Divide and Conquer MaxCut QAOA and Quantum Finance, demonstrate how quantum computing can be used outright to solve optimization problems. This notebook demonstrates a slightly different approach. Rather than considering QPUs as the tool to produce the final answer, it demonstrates how quantum can be used to enhance the effectiveness of leading classical methods.

The benefits of such an approach were highlighted in Scaling advantage with quantum-enhanced memetic tabu search for LABS. This notebook, co-created with the authors of the paper, will allow you to explore the findings of their research and write your own CUDA-Q code that builds a representative quantum-enhanced workflow for solving the LABS problem. Moreover, it will introduce advancements in counteradiabatic optimization techniques on which reduce the quantum resources required to run on a QPU.

Prerequisites: This lab assumes you have a basic knowledge of quantum computing, including operators, gates, etc. For a refresher on some of these topics, explore the Quick start to Quantum series.

In this lab you will:

1. Understand the LABS problem and its relation ot radar and communication applications.
1. Solve LABS classically with memetic tabu search and learn about the limitations of such methods.
1. Code a couteradiabatic algorithm using CUDA-Q to produce approximate solutions to the LABS problem.
1. Use the CUDA-Q results to seed your tabu search and understand the potential benefits of this approach.

Terminology you will use:

Low autocorrelation of binary sequences (LABS)
counteradiabatic optimization
memetic-tabu search

CUDA-Q Syntax you will use:

cudaq.sample()
@cudaq.kernel
ry(), rx(), rz(), x(), h()
x.ctrl()

The LABS problem and applications

The Low Autocorrelation Binary Sequences (LABS) problem is fundamental to many applications, but originated with applications to radar.

Consider a radar that monitors airport traffic. The radar signal sent to detect incoming planes must have as much range as possible to ensure safe approaches are planned well in advance. The range of a radar signal can be increased by sending a longer pulse. However, in order to differentiate between multiple objects, pulses need to be short to provide high resolution. So, how do you handle situations where you need both?

One solution is a technique called pulse compression. The idea is to send a long signal, but vary the phase at regular intervals such that the resolution is increased. Generally, the initial signal will encode a binary sequence of phase shifts, where each interval corresponds to a signal with a 0 or 180 degree phase shift.

The tricky part is selecting an optimal encoding sequence. When the signal returns, it is fed into a matched filter with the hope that a singular sharp peak will appear, indicating clear detection. The autocorrelation of the original signal, or how similar the signal is to itself, determines if a single peak or a messier signal with sidelobes will be detected. A signal should have high autocorrelation when overlayed on top of itself, but low autocorrelation when shifted with a lag.

Consider the image below. The signal on the left has a crisp single peak while the single on the right produces many undesirable sidelobes which can inhibit clear detection.

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.

E(s) = \sum_{k=1}^{N-1} C_k^2

Where $C_k$ is defined as.

C_k= \sum_{i=1}^{N-k} s_is_{i+k}

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. See if you can select a very good and very poor sequence and see the difference for the case of $N=7$ .

Classical Solution of the LABS problem

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.

Exercise 1:

Using the widget above, try to find some of the symmetries for the LABS problem. That is, for a fixed bitstring length, can you find patterns to produce the same energy with different pulse patterns.

The best known performance for a classical optimization technique is Memetic Tabu search (MTS) which exhibits a scaling of $O(1.34^N)$ . The MTS algorithm is depicted below. It begins with a randomly selected population of bitstrings and finds the best solution from them. Then, a child is selected by sampling directly from or combining multiple bitstrings from the population. The child is mutated with probability $p_{mutate}$ and then input to a tabu search, which performs a modified greedy local search starting from the child bitstring. If the result is better than the best in the population, it is updated as the new leader and randomly replaces a bitstring in the population.

Such an approach is fast, parallelizable, and allows for exploration with improved searching of the solution landscape.

Exercise 2:

Before exploring any quantum approach, get a sense for how MTS works by coding it yourself based on the figure above. Complete the TODO's below. Note that two functions are provided for you. `labs_utils` will perform the tabu search given a provided bitstring. `labs_plotting` will plot the energies of your final population and will optionally take as second data set when we later compare to a quantum enhanced result. The image below also defines the combine and mutate procedures you will need from the paper.

In [1]:

import sys
import os
sys.path.append('../auxiliary_files')
import numpy as np
import labs_utils as utils
import labs_plotting as plotting


def initialize_population(N, pop_size, user_population=None):
    """
    Initializes the population.
    If user_population is provided, it is used directly.
    Otherwise, a random population is generated.
    
    Args:
        N (int): Sequence length.
        pop_size (int): Number of individuals.
        user_population (np.array, optional): A specific population to start with.
        
    Returns:
        np.array: The initialized population.
    """
    if user_population is not None:
        return user_population
    else:
        #TODO - build initial random population


def combine(parent1, parent2):
    """
    Combination based on paper's algorithm 3
    The child takes the first part from p1 and the second part from p2.

    Args:
        parent_1 (np.array): sequence 1
        parent_2 (np.array): sequence 2
        
    Returns:
        np.array: Combined child of parent 1 and 2

    """
    #TODO Start - Write combine function

    #TODO End

def mutate(sequence, pmut):
    """
    Paper's Mutate function: Element-wise bit-flip with probability pmut.

    Args:
        sequence (np.array): input child sequence
        pmut (float): probability od mutation
        
    Returns:
        np.array: mutated string

    
    """
    #TODO Start -Write mutate function

    #TODO End

def run_memetic_search(initial_population, rounds_R):
    """
    Runs the Memetic Loop for R rounds.

    Args:
        rounds_R (int): number of memetic search rounds
        initial_population (np.array): array of sequences for the population bitstrings
        
    Returns:
        np.array: population energies
        float: best energy from population
    """

    #TODO Start - Write Memetic search algorithm
   
    #TODO END





# --- Configuration ---
N = 10         # Sequence Length
POP_SIZE = 100   # Population size
ROUNDS_R = 10   # How many memetic iterations (child generations) to run

print(f"--- LABS Memetic Search (N={N}, Rounds={ROUNDS_R}) ---\n")

print("Running with Random Initialization...")
pop_rand = initialize_population(N, POP_SIZE)

final_energies_rand, best_rand = run_memetic_search(pop_rand, ROUNDS_R)
print(f"    Best found: {best_rand}")


print("\nPlotting results...")
plotting.plot_energy_distributions(
    results_1=final_energies_rand, 
    labels=("Random Init + Memetic", "Quantum Enhanced + Memetic"),
    n_val=N
)

Out[1]:

--- LABS Memetic Search (N=10, Rounds=10) ---

Running with Random Initialization...
  > Starting Memetic Loop for 10 rounds...
  > Initial Best Energy: 13
    Best found: 13

Plotting results...

Building a Quantum Enhanced Workflow

Despite the effectiveness of MTS, it still exhibits exponential scaling $O(1.34^N)$ behavior and becomes intractable for large $N$ . Quantum computing provides a potential alternative method for solving the LABS problem because the properties of entanglement, interference, and superpositon may allow for a better global search. Recent demonstrations have even produced evidence that the quantum approximate optimization algorithm (QAOA) can be used to reduce the scaling of the LABS problem to $O(1.21^N)$ for $N$ between 28 and 40 with quantum minimum finding.

However, current quantum hardware limitations restrict solution to problems of greater than about $N=20$ , meaning that it will be some time before quantum approaches can outperform the classical state of the art. It should also be noted that standard QAOA can struggle with LABS and require many layers to converge the parameters if other tricks are not employed.

The authors of Scaling advantage with quantum-enhanced memetic tabu search for LABS cleverly explored an alternate path that combines quantum and classical approaches and might be able to provide a more near-term benefit. Instead of expecting the quantum computer to solve the problem entirely, they asked how a quantum approach might enhance MTS.

The basic idea is that a quantum optimization routine could run first and the resulting state be sampled to produce a better population for MTS. Many such heuristics for defining the initial population are possible, but the rest of this notebook will explore their methodology, help you to build the workflow yourself, and allow you to analyze the benefits of their approach.

The first step of quantum enhanced MTS (QE-MTS) is to prepare a circuit with CUDA-Q that approximates the ground state of the Hamiltonian corresponding to the LABS problem. You could do this with any optimization algorithm such as QAOA or using an adiabatic approach. (See the Quantum Portfolio Optimization CUDA-Q Academic lab for a detailed comparison of these two common approaches.)

The authors of this work opted for an adiabatic approach (More on why later). Recall that the goal of an adiabatic optimization is to begin with a Hamiltonian that has an easily prepared ground state ( $H_i$ ). Then, the adiabatic Hamiltonian $H_{ad}$ can be constructed as $H_{ad}(\lambda) = (1-\lambda)H_i +\lambda H_f$ , where $\lambda$ is a function of time and $H_f$ is the Hamiltonian representing a qubit encoding of the LABS problem.

H_f = 2 \sum_{i=1}^{N-2} \sigma_i^z \sum_{k=1}^{\lfloor \frac{N-i}{2} \rfloor} \sigma_{i+k}^z + 4 \sum_{i=1}^{N-3} \sigma_i^z \sum_{t=1}^{\lfloor \frac{N-i-1}{2} \rfloor} \sum_{k=t+1}^{N-i-t} \sigma_{i+t}^z \sigma_{i+k}^z \sigma_{i+k+t}^z

The authors also selected $H_i = \sum_i h^x_i \sigma^x_i$ which has an easily prepared ground state of $\ket{+}^{\otimes N}$ .

The challenge for implementing the optimization procedure becomes selection of an operator that will quickly and accurately evolve to the ground state of $H_f$ . One approach is to use a so-called auxiliary countradiabatic (CD) term $H_{CD}$ , which corrects diabatic transitions that jump out of the ground state during the evolution. The figure below demonstrates the benefit of using a CD correction.

An operator called the adiabatic gauge potential $A_{\lambda}$ is the ideal choice for the CD term as it suppresses all possible diabatic transitions, resulting in the following total system to evolve.

H(\lambda) = H_{ad}(\lambda) + \lambda H_{CD} (\lambda)

$A(\lambda)$ is derrived from $H_{ad}(\lambda)$ (see paper for details) as it contains information about underlying physics of the problem.

There is a problem though. The $A(\lambda)$ term cannot be efficiently expressed exactly and needs to be approximated. It also turns out that in the so called impulse regime, where the adiabatic evolution is very fast, $H_{cd} (\lambda)$ dominates $H_{ad}(\lambda)$ , meaning that the final implementation corresponds to the operator $H(\lambda) = H^1_{cd}(\lambda)$ where $H^1_{cd}(\lambda)$ is a first order approximation of $A(\lambda)$ (see equation 7 in the paper).

A final step is to use Trotterization to define the quantum circuit to apply $e^{-\theta (t) i H_{cd}}$ . The details for this derivation are shown in the appendix of the paper. and result from equation B3 is shown below.

\begin{equation} \begin{aligned} 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] \\ & \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) \\ & \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) \\ & \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) \\ & \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) \end{aligned} \end{equation}

It turns out that this implementation is more efficient than QAOA in terms of gate count. The authors calculated that for $N=67$ , QAOA would require 1.4 million entangling gates while the CD approach derived here requires only 236 thousand entangling gates.

Exercise 3:

At first glance, this equation might looks quite complicated. However, observe the structure and note two "blocks" of terms. Can you spot them?

They are 2 qubit terms that look like $R_{YZ}(\theta)$ or $R_{ZY}(\theta)$ .

As well as 4 qubit terms that look like $R_{YZZZ}(\theta)$ , $R_{ZYZZ}(\theta)$ , $R_{ZZYZ}(\theta)$ , or $R_{ZZZY}(\theta)$ .

Thankfully the authors derive a pair of circuit implementations for the two and four qubit terms respectively, shown in the figures below.

Using CUDA-Q, write a kernel for each which will be used later to construct the full implementation.

Hint: Remember that the adjoint of a rotation gate is the same as rotating in the opposite direction.
Hint: You may also want to define a CUDA-Q kernel for an R $_{ZZ}$ gate.
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.

In [2]:

import cudaq 
@cudaq.kernel
def rzz_gate(angle: float, q0: cudaq.qubit, q1: cudaq.qubit):

    #TODO - write kernel for RZZ gate



@cudaq.kernel
def two_qubit_block(angle: float, q0: cudaq.qubit, q1: cudaq.qubit):
    
    # TODO - Write kernel for two qubit blocks


@cudaq.kernel
def four_qubit_block(angle: float, q0: cudaq.qubit, q1: cudaq.qubit, q2: cudaq.qubit, q3: cudaq.qubit):

     # TODO - Write kernel for four qubit blocks

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.

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)$ .

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)$ .

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

\alpha(t) = \frac{-\Gamma_1(t)}{\Gamma_2(t)}

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.

Exercise 4:

The `compute_theta` function below 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 list 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.

\begin{equation} \begin{aligned} 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] \\ & \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) \\ & \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) \\ & \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) \\ & \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) \end{aligned} \end{equation}

In [3]:

from math import floor

#TODO - complete the loops below to save G2 and G4 indicies
def get_interactions(N):
    """
    Generates the interaction sets G2 and G4 based on the loop limits in Eq. 15.
    Returns standard 0-based indices as lists of lists of ints.
    
    Args:
        N (int): Sequence length.
        
    Returns:
        G2: List of lists containing two body term indices
        G4: List of lists containing four body term indices
    """
    #TODO - complete the loops below to save G2 and G4 indicies
    G2 = []
    
 

    G4 = []


#TODO END
                
    return G2, G4

Exercise 5:

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.

Use cudaq.sample to extract the samples we want. You will need to post-process these to the appropriate format for input to the MTS population.

In [4]:



@cudaq.kernel
def trotterized_circuit(N: int, G2: list[list[int]], G4: list[list[int]], steps: int, dt: float, T: float, thetas: list[float]):
    
    #TODO - Write kenrel to apply full circuit


  

T=1
n_steps = 1
dt = T / n_steps
N = 20
G2, G4 = get_interactions(N)

thetas =[]

for step in range(1, n_steps + 1):
    t = step * dt
    theta_val = utils.compute_theta(t, dt, T, N)
    thetas.append(theta_val)

print(cudaq.sample(trotterized_circuit, N, G2, G4, n_steps, dt, T, thetas, shots_count=10))

Out[4]:

{ 00011100001101010001:1 00101001101010010000:1 01001101100000110000:1 01100110000101011011:1 01110001010111001010:1 10001010100101000011:1 10010100111111010010:1 10011010000000000001:1 10101011010001111000:1 10110011000001111000:1 }

Generating Quantum Enhanced Results

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.

Exercise 6:

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 have a GPU available, try to run a larger problem size of at least $N=20$ . Feel free to explore other problem sizes too. You can still run the cell below with a CPU, but change the sequence length to something small like $N=5$ .

In [5]:

T=1
n_steps = 1
dt = T / n_steps
N = 20  # Sequence Length
POP_SIZE = 100   # Population size
ROUNDS_R = 10   # How many memetic iterations (child generations) to run
G2, G4 = get_interactions(N)

thetas =[]

for step in range(1, n_steps + 1):
    t = step * dt
    theta_val = utils.compute_theta(t, dt, T, N)
    thetas.append(theta_val)


samples = cudaq.sample(trotterized_circuit, N, G2, G4, n_steps, dt, T, thetas, shots_count=POP_SIZE)



samples_list = []
for bitstring, count in samples.items():
    bits = [int(b) for b in bitstring]
    samples_list.extend([bits] * count)

samples = np.array(samples_list)


print(f"--- LABS Memetic Search (N={N}, Rounds={ROUNDS_R}) ---\n")

print("Running with Random Initialization...")
pop_rand = initialize_population(N, POP_SIZE)

final_energies_rand, best_rand = run_memetic_search(pop_rand, ROUNDS_R)
print(f"    Best found Random Population: {best_rand}")

final_energies_quantum, best_quantum = run_memetic_search(samples, ROUNDS_R)
print(f"    Best found Quantum Population: {best_quantum}")

print("\nPlotting results...")
plotting.plot_energy_distributions(
    results_1=final_energies_rand, 
    results_2=final_energies_quantum,
    labels=("Random Init + Memetic", "Quantum Enhanced + Memetic"),
    n_val=N
)

Out[5]:

--- LABS Memetic Search (N=20, Rounds=10) ---

Running with Random Initialization...
  > Starting Memetic Loop for 10 rounds...
  > Initial Best Energy: 62
    Best found Random Population: 38
  > Starting Memetic Loop for 10 rounds...
  > Initial Best Energy: 6
    Best found Quantum Population: 5

Plotting results...

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.

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.

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.

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.

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.

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.

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)$ ).

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.

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.

In [ ]: