# The Maximum Independent Set Problem

## Background

In graph theory, an independent set is a set of vertices in a graph such that no two of which are connected by an edge. The problem of finding Maximum Independent Sets (MISs) is NP-hard, i.e., it is unlikely to be solved in a time polynomial to the problem size. Interestingly, there is a natural connection between the independent set constraint, and the Rydberg Blockade phenomenon in neutral-atom quantum computing using Rydberg states. More specifically, the Rydberg blockade implies that two atoms cannot be both excited to the Rydberg state $\Ket{r}$ if they are close to each other, whereas the independent set constraint means two vertices cannot be both in the independent set when they are connected by an edge. Thus, one can consider atoms in the Rydberg state as vertices in an independent set. See the proposal in H. Pichler et al. (10.48550/arXiv.1808.10816) for more details.

In particular, one can use the ground state of the Rydberg Hamiltonian to encode the Maximum Independent Set problem, which is to find the largest independent set of a given graph. For a particular subclass of geometric graphs, the so-called unit disk graphs, the Rydberg Hamiltonian can encode the solution without any overhead in the number of qubits. In fact, an experimental demonstration of quantum optimization has been realized in solving the Maximum Independent Set problem on up to 289 qubits in S. Ebadi et al. (10.48550/arXiv.2202.09372).

In this tutorial, we show how to solve the MIS problem using Bloqade. We focus on a particular subclass of unit disk graphs defined as Diagonal-connected Unit-disk Grid Graphs (DUGG). This is the class of graphs studied in the demonstration experiment S. Ebadi et al. (10.48550/arXiv.2202.09372). Although these graphs have highly constrained topologies, finding its MISs is still NP-hard. Here, we show how to use variational quantum algorithms with Rydberg Hamiltonians to solve the MIS problem on these graphs. The tutorial here strongly resembles the setup in S. Ebadi et al. (10.48550/arXiv.2202.09372), with the limitation that Bloqade can only simulate a much smaller problem.

For more details on the functionalities supported by Bloqade in studying independent set problems, please refer to the Maximum Independent Set page in the Manual.

Let's start by importing the required libraries:

```
using Graphs
using Bloqade
using Random
using GenericTensorNetworks
using Optim
using PythonCall
plt = pyimport("matplotlib.pyplot");
```

## Setting Up the Problem

To begin, we create a $4 \times 4$ DUGG with 0.8 filling, by using the `random_dropout`

function. Here, we choose the lattice constant $a$ to be 4.5 μm:

```
Random.seed!(2)
atoms = generate_sites(SquareLattice(), 4, 4; scale = 4.5) |> random_dropout(0.2)
```

Next, we set the blockade radius to be 7.5 μm, corresponding to a case where nearest neighbors and next-nearest neighbors (diagonal from the initial atom) are within the blockade radius. As we discussed in the Rydberg Blockade tutorial, only one Rydberg excitation is allowed within the blockade radius. To better illustrate this constraint, we plot the interactions of Rydberg atoms as a DUGG, where each edge corresponds to the blockade constraint given by the blockade radius:

`Bloqade.plot(atoms, blockade_radius = 7.5)`

Our goal is to find a maximum independent set of such a graph.

For comparison, we first use classical algorithms to calculate the MIS size here using the graph utilities in Bloqade, so that one can compare this exact result with the quantum algorithms. The exact MIS size and its degeneracy can be solved with the generic tensor network algorithm (J. Liu et al. (10.48550/arXiv.2205.03718)) in the package `GenericTensorNetworks`

:

```
graph = BloqadeMIS.unit_disk_graph(atoms, 7.5)
mis_size_and_counting = GenericTensorNetworks.solve(IndependentSet(graph), CountingMax())[]
```

`(4.0, 26.0)ₜ`

The `solve`

function takes a graph instance and a solution space property as inputs, where the graph instance is generated by the `unit_disk_graph`

function in the Bloqade submodule `BloqadeMIS`

, and the solution space property is to count the number of maximum independent sets here. For this specific DUGG, we see that the MIS size is 4, and the number of possible Maximum Independent Sets is 26. In the following, we are going to show how to solve the independent set problem with both quantum adiabatic and variational algorithms.

## The Adiabatic Approach

Here, we generalize the quantum adiabatic algorithm used in the Adiabatic Evolution tutorial to prepare ground states of the Rydberg Hamiltonian for this disordered lattice. We first construct the adiabatic pulse sequences for the Rabi frequency $\Omega$ and the detuning $\Delta$:

```
T_max = 0.6
Ω_max = 2π * 4
Ω = piecewise_linear(clocks = [0.0, 0.1, 0.5, T_max], values = [0.0, Ω_max, Ω_max, 0])
Δ_start = -2π * 13
Δ_end = 2π * 11
Δ = piecewise_linear(clocks = [0.0, 0.1, 0.5, T_max], values = [Δ_start, Δ_start, Δ_end, Δ_end])
fig, (ax1, ax2) = plt.subplots(ncols = 2, figsize = (12, 4))
Bloqade.plot!(ax1, Ω)
ax1.set_ylabel("Ω/2π (MHz)")
Bloqade.plot!(ax2, Δ)
ax2.set_ylabel("Δ/2π (MHz)")
fig
```

Here, the total time is fixed to `T_max`

, the adiabatic evolution path is specified by the `piecewise_linear`

function. The Rydberg blockade radius can be computed with

\[C_6 / R_b^6 \sim \sqrt{\Delta^2 + \Omega^2}\]

The default for $C_6=2π * 862690 \text{ MHz μm}^6$. For encoding the corresponding MIS problem at $\Omega = 0$, the detuning can be set around $2\pi \times 11$ MHz for a blockade radius of $7.5$ µm (see the parameters in S. Ebadi et al. (10.48550/arXiv.2202.09372)).

Next, we create the time-dependent Hamiltonian and simulate its time evolution by using the `SchrodingerProblem`

solver:

```
hamiltonian = rydberg_h(atoms; Ω = Ω, Δ = Δ)
prob = SchrodingerProblem(zero_state(nqubits(hamiltonian)), T_max, hamiltonian)
emulate!(prob)
```

```
SchrodingerProblem:
register info:
type: YaoArrayRegister.ArrayReg{2, ComplexF64, Matrix{ComplexF64}}
storage size: 8 bytes
time span (μs): (0.0, 0.6)
equation:
storage size: 1.688 MiB
expression:
nqubits: 13
+
├─ [+] ∑ 2π ⋅ 8.627e5.0/|x_i-x_j|^6 n_i n_j
├─ [+] Ω(t) ⋅ ∑ σ^x_i
└─ [-] Δ(t) ⋅ ∑ n_i
algorithm: DP8(; stage_limiter! = trivial_limiter!, step_limiter! = trivial_limiter!, thread = static(false),)
options:
save_everystep: false
save_start: false
save_on: false
dense: false
reltol: 1.0e-10
abstol: 1.0e-10
```

Finally, we can plot the most probable bitstrings by using the `bitstring_hist`

function on the resulting register (quantum state):

`bitstring_hist(prob.reg; nlargest = 20)`

We can see that some of the most probable configurations indeed have an independent set size 4 by counting the number of ones in the bitstring. The correctness of the output can be verified by comparing it to the classical solution:

```
best_bit_strings = most_probable(prob.reg, 2)
all_optimal_configs = GenericTensorNetworks.solve(IndependentSet(graph), ConfigsMax())[]
@assert all(bs -> GenericTensorNetworks.StaticBitVector([bs...]) ∈ all_optimal_configs.c, best_bit_strings)
```

We can also visualize these atoms and check them visually:

`Bloqade.plot(atoms, blockade_radius = 7.5; colors = [iszero(b) ? "white" : "red" for b ∈ best_bit_strings[1]])`

`Bloqade.plot(atoms, blockade_radius = 7.5; colors = [iszero(b) ? "white" : "red" for b ∈ best_bit_strings[2]])`

However, there are still some configurations that violate the blockade constraint, because the blockade interaction is not an ideal unit disk constraint (e.g. some bitstrings have a size 5). One can check whether the independence constraint is satisfied or not with the `BloqadeMIS.is_independent_set`

function:

```
best5_bit_strings = most_probable(prob.reg, 3)
BloqadeMIS.is_independent_set.(best5_bit_strings, Ref(graph))
```

```
3-element BitVector:
1
1
0
```

This can happen when the Rydberg interaction is not very strong at the unit disk radius. See the Rydberg Blockade page for more detailed explanations on Rydberg blockade and its relation with the unit disk radius. One can perform some postprocessing by reducing the violated configurations to indendendent set configurations using the `mis_postprocessing`

function. Please also refer to the paper S. Ebadi et al. (10.48550/arXiv.2202.09372)). for more detailed discussion on the postprocessing procedure.

```
fixed = mis_postprocessing(best5_bit_strings[3], graph)
BloqadeMIS.is_independent_set(fixed, graph)
```

`true`

## QAOA with Piecewise Constant Pulses

The QAOA (Quantum Approximate Optimization Algorithm) algorithm (E. Farhi, J. Goldstone, S. Gutmann (10.48550/arXiv.1411.4028)) is a hybrid quantum-classical algorithm. The classical part of the algorithm is an optimizer, which can be either a gradient-based or non-gradient-based one. For our specific problem, the corresponding quantum part is a neutral-atom quantum computer first evolving under the Rydberg Hamiltonian with parameterized pulse sequences and then being measured in the computational basis.

The standard definition of QAOA involves applying the problem (expressed as a cost function) Hamiltonian $C$ and the transverse field Hamiltonian $B$ alternately. Let $G=(V,E)$ be a graph. The cost Hamiltonian for an MIS problem can be defined as

\[C(G) = -\sum_{j\in V}^{n} w_j n_j^z + \infty \sum_{\langle j,k\rangle \in E}n_j^z n_k^z\]

where the first summation is proportional to the size of the independent set, and the second term enforces the independence constraints.

In a Rydberg Hamiltonian, the first term corresponds to the detuning $w_i = \Delta$. The second term contains an $\infty$, which corresponds to the Rydberg blockade term with its strength described as $V_{jk} = C_6/|\mathbf{x}_j - \mathbf{x}_k|^6$. As we can see, the Rydberg interaction is not a perfect independence constraint due to with finite blockade interaction and unwanted long-range interaction. Thus, postprocessing might be required using neutral-atom quantum computers to solve the MIS problem.

The transverse-field Hamiltonian corresponds to the Rabi term in the Rydberg Hamiltonian,

\[B = \sum_{j=1}^{n}\sigma_j^x + \infty \sum_{\langle j,k\rangle \in E}n_j^z n_k^z.\]

Note that the Rybderg interaction term is always on, in contrast to the standard QAOA protocol. For the convenience of the simulation, we use the `expect`

function to get the averaged measurement outputs. On the actual quantum hardware, the `expect`

should be replaced by measuring in the computational basis and obtaining the averaged number of Rydberg excitations as the loss function (also called objective function or cost function). One can then either use non-gradient-based optimizers to perform the optimization or finite-difference methods to obtain gradients of parameters.

Let us first set up a non-optimized pulse sequence for QAOA with step $p=3$:

```
durations = fill(0.1, 6)
clocks = [0, cumsum(durations)...]
Ω2 = piecewise_constant(; clocks = clocks, values = repeat([Ω_max, 0.0], 3))
Δ2 = piecewise_constant(; clocks = clocks, values = repeat([0.0, Δ_end], 3))
fig, (ax1, ax2) = plt.subplots(ncols = 2, figsize = (12, 4))
Bloqade.plot!(ax1, Ω2)
ax1.set_ylabel("Ω/2π (MHz)")
Bloqade.plot!(ax2, Δ2)
ax2.set_ylabel("Δ/2π (MHz)")
fig
```

The `piecewise_constant`

pulses can be more accurately simulated with the `KrylovEvolution`

solver. This time, we simulate the dynamics in the subspace generated by the `blockade_subspace`

function, so that we do not need postprocessing anymore.

```
hamiltonian2 = rydberg_h(atoms; Ω = Ω2, Δ = Δ2)
nsites = length(atoms)
subspace = blockade_subspace(atoms, 7.5) # we run our simulation within the blockade subspace
prob2 = KrylovEvolution(zero_state(subspace), clocks, hamiltonian2)
emulate!(prob2);
```

We defined the loss function as the negative of the mean MIS size, which corresponds to the expectation value of the `SumOfN`

operator. Thus, we can calculate the average loss function after the time evolution:

```
loss_MIS(reg) = -rydberg_density_sum(reg)
loss_MIS(prob2.reg)
```

`-2.5628869128026794`

The output shows the negative mean independent set size. This is because we have flipped its sign since most optimizers are set to minimize the loss function. This loss is equivalent to the `rydberg_density_sum`

loss function in the `BloqadeMIS`

module. Alternative loss functions include the `gibbs_loss`

and the `independent_set_probabilities`

.

Here, the loss produced by the pulse sequence does not look very good. We can throw it into an optimizer and see if a classical optimizer can help. First, let us wrap up the above code into a loss function:

```
function loss_piecewise_constant(atoms::AtomList, x::AbstractVector{T}) where {T}
@assert length(x) % 2 == 0
Ω_max = 4 * 2π
Δ_end = 11 * 2π
p = length(x) ÷ 2
# detuning and rabi terms
durations = abs.(x) # the durations of each layer of the QAOA pulse take the optimizing vector x as their input
clocks = [0, cumsum(durations)...]
Ωs = piecewise_constant(; clocks = clocks, values = repeat(T[Ω_max, 0.0], p))
Δs = piecewise_constant(; clocks = clocks, values = repeat(T[0.0, Δ_end], p))
hamiltonian = rydberg_h(atoms; Ω = Ωs, Δ = Δs)
subspace = blockade_subspace(atoms, 7.5) # we run our simulation within the blockade subspace
prob = KrylovEvolution(zero_state(Complex{T}, subspace), clocks, hamiltonian)
emulate!(prob)
return -rydberg_density_sum(prob.reg), prob.reg
end
```

`loss_piecewise_constant (generic function with 1 method)`

Running the simulation in the subspace does not violate the independence constraints. If one uses fullspace simulation or runs it on the quantum computer, one may need to post-process the measured bitstrings to a get a correct measure of the loss if we do not set the blockade constraint. Related APIs include `is_independent_set`

, `num_mis_violation`

, and `mis_postprocessing`

.

Let us check the loss function again using the initial point above:

```
x0 = durations
rydberg_density, reg1 = loss_piecewise_constant(atoms, x0)
rydberg_density
```

`-2.5628869128026794`

The most probable bitstrings are:

`bitstring_hist(reg1; nlargest = 20)`

We see that without optimization, many of these bitstrings are not the MIS solutions.

Let us now use the non-gradient-based optimizer `NelderMead`

in the `Optim`

package to optimize the loss function:

```
optresult = Optim.optimize(x -> loss_piecewise_constant(atoms, x)[1], x0)
rydberg_density_final, reg1_final = loss_piecewise_constant(atoms, optresult.minimizer)
rydberg_density_final
```

`-3.0965910260012097`

We see that the loss is indeed improved, but still not very good. This is likely because the optimization is trapped in a local minimum:

`bitstring_hist(reg1_final; nlargest = 20)`

This example shows that the performance of the algorithm very much depends on the parametrization of the pulse sequences, the initialization of the variational parameters, and the classical optimizers. See S. Ebadi et al. (10.48550/arXiv.2202.09372) for more in-depth comparison of different pulse parameterizations and tips on how to improve the performance.

In the example below, we show a better pulse parameterization using smoothened piecewise linear waveforms.

## Smoothened Piecewise Linear Pulses

A smoothened piecewise linear waveform can be created by applying a Gaussian filter on a waveform created by the `piecewise_linear`

function. For example:

```
pulse_piecewise_linear = piecewise_linear(clocks = [0.0, 0.05, 0.1, 0.5, 0.55, T_max], values = [0, 0, 0.4, 0.4, 0, 0]);
pulse_smooth = smooth(pulse_piecewise_linear; kernel_radius = 0.02);
fig, ax = plt.subplots()
Bloqade.plot!(ax, pulse_piecewise_linear)
Bloqade.plot!(ax, pulse_smooth)
ax.set_ylabel("strength")
ax.legend(["piecewise linear", "smoothened piecewise linear"], loc = "lower right")
fig
```

Here, the function `smooth`

takes a `kernel_radius`

keyword parameter as the Gaussian kernel parameter. With the new waveforms, we can define the loss function as follows:

```
function loss_piecewise_linear(atoms::AtomList, x::AbstractVector{T}) where {T}
@assert length(x) == 3
Ω_max = 4 * 2π
Δ_start = -13 * 2π
Δ_end = 11 * 2π
Δ0 = 11 * 2π
T_max = 0.6
# the strength of the detunings at each step takes the optimizing x as their input
Δs = smooth(
piecewise_linear(
clocks = T[0.0, 0.05, 0.2, 0.3, 0.4, 0.55, T_max],
values = T[Δ_start, Δ_start, Δ0*x[1], Δ0*x[2], Δ0*x[3], Δ_end, Δ_end],
);
kernel_radius = 0.02,
)
Ωs = smooth(
piecewise_linear(clocks = T[0.0, 0.05, 0.1, 0.5, 0.55, T_max], values = T[0, 0, Ω_max, Ω_max, 0, 0]);
kernel_radius = 0.02,
)
hamiltonian = rydberg_h(atoms; Ω = Ωs, Δ = Δs)
subspace = blockade_subspace(atoms, 7.5)
prob = SchrodingerProblem(zero_state(Complex{T}, subspace), T_max, hamiltonian)
emulate!(prob)
return -rydberg_density_sum(prob.reg), prob.reg, Δs
end
x0 = [0.1, 0.8, 0.8]; # initial point for the optimization
```

Let us check the loss function with smoothened waveform with the initial point:

```
Δ_start = -13 * 2π
Δ_end = 11 * 2π
Δ0 = 11 * 2π
T_max = 0.6
Δ_initial = piecewise_linear(
clocks = [0.0, 0.05, 0.2, 0.3, 0.4, 0.55, T_max],
values = [Δ_start, Δ_start, Δ0 * x0[1], Δ0 * x0[2], Δ0 * x0[3], Δ_end, Δ_end],
)
rydberg_density, reg2, Δ_initial_smooth = loss_piecewise_linear(atoms, x0)
rydberg_density
```

`-3.890385062164213`

and plot the smoothened waveform:

```
fig, ax = plt.subplots()
Bloqade.plot!(ax, Δ_initial)
Bloqade.plot!(ax, Δ_initial_smooth)
ax.set_ylabel("Δ/2π (MHz)")
ax.legend(["piecewise linear", "smoothened piecewise linear"], loc = "lower right")
fig
```

Let's plot the distribution:

`bitstring_hist(reg2; nlargest = 20)`

The performance of the algorithm is quite good. Again, let us use the `NelderMead`

optimizer to optimize the loss function:

```
optresult = Optim.optimize(x -> loss_piecewise_linear(atoms, x)[1], x0)
rydberg_density_final, reg_final, Δ_final = loss_piecewise_linear(atoms, optresult.minimizer)
rydberg_density_final
```

`-3.974827326383817`

One can see the mean MIS size can be further improved to a value close to the size of the MIS, which means there is a substantial probability for measuring an MIS state.

`bitstring_hist(reg_final; nlargest = 20)`

We can also plot out the final optimized waveform for Δ and compare it with the initial waveform:

```
fig, ax = plt.subplots()
Bloqade.plot!(ax, Δ_initial_smooth)
Bloqade.plot!(ax, Δ_final)
ax.set_ylabel("Δ/2π (MHz)")
ax.legend(["initial", "optimized"], loc = "lower right")
fig
```

## Applications in VLSI Chip Manufacturing

### Background

The Maximum Independent Set (MIS) tutorial above described how to use Bloqade to computed the MIS on arbitrary unit disk graphs. In this tutorial, we will use the same algorithm but also show how to map a problem with real-world applications to a problem that can be solved with Bloqade's features.

Configurability is an approach to chip manufacturing in VLSI design. Ordinarily, copies of the entire chip and fabricated and on a wafer, the wafer is diced into chips using scribing corridors, and then faulty chips are marked and discarded. However, since the entire chip is being manufactured at once, this is wasteful because even if only a small portion of the chip is defective, it must be entirely discarded. According to F. Berman et al., "To use configurability, this process is changed somewhat. Special circuitry is added in the scribing corridors that enables adjacent chips to be connected. Then, after the wafer is tested, adjacent functional chips are connected into a group and the entire group is used as a single enlarged chip." In this way, even if a smaller chips becomes defective, the adjacent chips can be rewired so that they can still be used in the larger chip.

A natural problem that arises from this setup is the *square packing* problem. D. Hochbaum and W. Maass summarizes an example problem: "64K RAM chips, some of which may be defective, are available on a rectilinear grid placed on a silicon wafer. 2 x 2 arrays of such nondefective chips could be wired together to produce 256K RAM chips. In order to maximize yield, we want to pack a maximal number of such 2 X 2 arrays into the array of working chips on a wafer."

To solve this problem, let us start with the necessary imports:

```
using Bloqade
using GenericTensorNetworks
using GenericTensorNetworks: unit_disk_graph
using PythonCall
using Random
using Optim
plt = pyimport("matplotlib.pyplot");
patches = pyimport("matplotlib.patches");
np = pyimport("numpy")
```

`Python: <module 'numpy' from '/home/runner/work/Bloqade.jl/Bloqade.jl/examples/1.blockade/.CondaPkg/env/lib/python3.10/site-packages/numpy/__init__.py'>`

### Setting up the Problem

We can visualize this problem as a grid of 64K RAM chips, where some of the grid cells are removed.

```
fig, ax = plt.subplots(figsize=(5, 5))
G = 5 # grid size (G x G)
defects = []
for i in range(1, G*G//15)
push!(defects, ((rand(Int)%G+10) % G, (rand(Int)%G + 10) % G))
end
ax.grid(visible=true)
plt.xticks(1:G)
plt.yticks(1:G)
grid = zeros((G, G))
function plot_defects(defects)
for (x, y) in defects
ax.add_patch(patches.Rectangle((x, y), 1, 1, color="red"))
grid[x+1, y+1] = 1
end
end
plot_defects(defects)
fig
```

Wherever there is a 2x2 grid of non-defective chips, these grid cells can be wired together. This potential grouping of chips can be visualized as a unit disk which overlaps with the four grid cells that are being wired together. For example, in the figure below, the four grid cells in the upper left corner could be wired together, so we would put a unit disk at (1, 9) with radius 1.

The code below shows the unit disk corresponding to every possible wiring. Overlapping unit disks represents two wirings which conflict with each other because they use the same chips. Our goal is to find the maximum number of non-conflicting wirings, i.e. the MIS.

```
circs::Vector{Tuple{Int64, Int64}} = []
Δ = [(0, 0) (0, 1) (1, 0) (1, 1)]
for x in 1:G-1
for y in 1:G-1
works = true
for (dx, dy) in Δ
if 1 ≤ x+dx ≤ G && 1 ≤ y + dy ≤ G
works &= grid[x+dx,y+dy]==0
end
end
if works
push!(circs, (x, y))
end
end
end
R_opt = sqrt(2sqrt(2))/2
for circ in circs
patch = patches.Circle(circ, R_opt, fill=false)
ax.add_patch(patch)
end
fig
```

We can now use Bloqade, transforming this problem (which is a type of Diagonal-connected Unit-disk Grid Graph (DUGG)) into a lattice of atoms. We want each atom to be positioned at the center of the unit disk and have a blockade radius of roughly 7.5 μm, reusing the methdology and values given in the "Setting Up the Problem" section of the tutorial above which enables atoms that exist diagonally from the initial atom to fall within the Blockade radius.

We will scale the position of the atoms so that we get the desired unit disk overlaps with a radius of 7.5 μm:

```
function scalePos(pos, scale)::Vector{Tuple{Float64, Float64}}
npos = []
for (x, y) in pos
push!(npos, (x*scale, y*scale))
end
return npos
end
R = 7.5
scale = 7.5/(2*R_opt)
atoms = AtomList(scalePos(circs, scale))
Bloqade.plot(atoms, blockade_radius = R)
```

We will first use classical algorithms to find the optimal solution to the problem:

```
unit_graph = unit_disk_graph(atoms, R)
configs = GenericTensorNetworks.solve(IndependentSet(unit_graph), ConfigsMax())[]
MIS_config = configs.c[1]
```

`01100000001001`

We can visualize the solution on the atom lattice:

`Bloqade.plot(atoms, blockade_radius = R; colors = [iszero(b) ? "white" : "red" for b in MIS_config])`

We can also visualize the solution on the chip grid:

```
function clear_circ()
for i in 1:5
for patch in ax.patches
try
patch.radius
patch.remove()
catch
print(patch)
end
end
end
end
clear_circ()
for (circ, used) in zip(circs, MIS_config)
if used == 1
patch = patches.Circle(circ, R_opt, fill=false)
ax.add_patch(patch)
end
end
fig
```

### The Adiabatic Approach

Like with the Maximum Independent Set tutorial above, we can use an adiabatic algorithm to solve this problem. We can use the same pulse sequence for $\Omega$ and $\Delta$:

```
T_max = 0.6
Ω_max = 2π * 4
Ω = piecewise_linear(clocks = [0.0, 0.1, 0.5, T_max], values = [0.0, Ω_max, Ω_max, 0])
Δ_start = -2π * 13
Δ_end = 2π * 11
Δ = piecewise_linear(clocks = [0.0, 0.1, 0.5, T_max], values = [Δ_start, Δ_start, Δ_end, Δ_end])
plot, (ax1, ax2) = plt.subplots(ncols = 2, figsize = (12, 4))
Bloqade.plot!(ax1, Ω)
ax1.set_ylabel("Ω/2π (MHz)")
Bloqade.plot!(ax2, Δ)
ax2.set_ylabel("Δ/2π (MHz)")
plot
```

And construct the Hamiltonian:

```
hamiltonian = rydberg_h(atoms; Ω = Ω, Δ = Δ, ϕ=Waveform(t->0.0, T_max))
prob = SchrodingerProblem(zero_state(nqubits(hamiltonian)), T_max, hamiltonian)
emulate!(prob)
```

```
SchrodingerProblem:
register info:
type: YaoArrayRegister.ArrayReg{2, ComplexF64, Matrix{ComplexF64}}
storage size: 8 bytes
time span (μs): (0.0, 0.6)
equation:
storage size: 5.500 MiB
expression:
nqubits: 14
+
├─ [+] ∑ 2π ⋅ 8.627e5.0/|x_i-x_j|^6 n_i n_j
├─ [+] Ω(t) ⋅∑ e^{ϕ(t) ⋅ im} |0⟩⟨1| + e^{-ϕ(t) ⋅ im} |1⟩⟨0|
└─ [-] Δ(t) ⋅ ∑ n_i
algorithm: DP8(; stage_limiter! = trivial_limiter!, step_limiter! = trivial_limiter!, thread = static(false),)
options:
save_everystep: false
save_start: false
save_on: false
dense: false
reltol: 1.0e-10
abstol: 1.0e-10
```

We can see that our algorithm found the Maximum Independent Set by summing the probabilities of being in any one of the MIS configurations:

```
function get_MIS_prob(reg::Union{ArrayReg, SubspaceArrayReg}, configs) # want to maxmimize this
prob = 0
x = [parse(Int, reverse(string(x)); base=2) for x in configs.c]
for (c, amp) in BloqadeMIS.ConfigAmplitude(reg)
if c in x
prob+=abs2(amp)
end
end
return prob
end
get_MIS_prob(prob.reg, configs)
```

`0.8652655503921888`

### Optimization

For even better results, we can optimize the pulse sequence for $\Omega$ and $\Delta$. We can parametrize the pulses by $t_{start}$ and $t_{end}$, which are the times spent in the initial and final state, respectively. We will also optimize the shape on the detuning pulse by modeling it as a piecewise linear function given by parameters that indicate what values it reaches at regular time intervals. We will optimize for the probabilitity that of getting one the MIS solutions without applying any classical postprocessing on the results of the stimulation.

```
function get_waves(params::Vector{Float64}) # returns the pulses as a function of certain parameters
t_start, t_end, scale... = params
t_start+=0.03
T_max = 0.6
t_end = T_max - t_end - 0.03
t_interval = (t_end - t_start) / (length(scale) + 1)
Ω_max = 4 * 2π
Δ_start = -13 * 2π
Δ_end = Δ0 = 11 * 2π
Δ_clock = [0.0, t_start]
Δ_val = [Δ_start, Δ_start]
for i in 1:length(scale)
push!(Δ_val, Δ0*scale[i])
push!(Δ_clock, t_start + i*t_interval)
end
push!(Δ_val, Δ_end); push!(Δ_val, Δ_end)
push!(Δ_clock, t_end); push!(Δ_clock, T_max)
Δ = piecewise_linear(clocks = Δ_clock, values = Δ_val);
Ω = piecewise_linear(clocks = [0.0, t_start, t_end, T_max], values = [0, Ω_max, Ω_max, 0]);
return Ω, Δ
end
```

`get_waves (generic function with 1 method)`

The example below plots the pulse sequence that would be generated by the parameters given in `x0`

:

```
function plot_waves(params::Vector{Float64})
Ω, Δ = get_waves(params)
graph, (ax1, ax2) = plt.subplots(ncols = 2, figsize = (12, 4))
Bloqade.plot!(ax1, Ω)
ax1.set_ylabel("Ω/2π (MHz)")
Bloqade.plot!(ax2, Δ)
ax2.set_ylabel("Δ/2π (MHz)")
return graph
end
x0 = [0.1, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8]
plot_waves(x0)
```

We can now run the optimization using the `Nelder-Mead`

algorithm

```
function loss(atoms::AtomList, x::Vector{Float64}, configs)
Ω, Δ = get_waves(x)
hamiltonian = rydberg_h(atoms; Ω = Ω, Δ = Δ)
subspace = blockade_subspace(atoms, 4.5)
prob = SchrodingerProblem(zero_state(subspace), T_max, hamiltonian)
emulate!(prob)
return -get_MIS_prob(prob.reg, configs), prob.reg, Δ
end
optresult = Optim.optimize(x -> loss(atoms, x, configs)[1], x0)
```

```
* Status: success
* Candidate solution
Final objective value: -9.676880e-01
* Found with
Algorithm: Nelder-Mead
* Convergence measures
√(Σ(yᵢ-ȳ)²)/n ≤ 1.0e-08
* Work counters
Seconds run: 85 (vs limit Inf)
Iterations: 463
f(x) calls: 735
```

We can look at the new pulse sequence, which has a much better `get_MIS_prob`

:

`plot_waves(optresult.minimizer)`

The new parameters give much better results when we rerun the adiabtic algorithm:

```
Ω, Δ = get_waves(optresult.minimizer)
hamiltonian = rydberg_h(atoms; Ω = Ω, Δ = Δ, ϕ=Waveform(t->0.0, T_max))
prob = SchrodingerProblem(zero_state(nqubits(hamiltonian)), T_max, hamiltonian)
emulate!(prob)
bitstring_hist(prob.reg; nlargest = 20)
```

*This page was generated using Literate.jl.*