Hardness of the maximum-independent-set problem on unit-disk graphs and prospects for quantum speedups

Rydberg atom arrays are among the leading contenders for the demonstration of quantum speedups. Motivated by recent experiments with up to 289 qubits [Ebadi et al., Science 376, 1209 (2022)], we study the maximum-independent-set problem on unit-disk graphs with a broader range of classical solvers beyond the scope of the original paper. We carry out extensive numerical studies and assess problem hardness, using both exact and heuristic algorithms. We find that quasiplanar instances with Union-Jack-like connectivity can be solved to optimality for up to thousands of nodes within minutes, with both custom and generic commercial solvers on commodity hardware, without any instance-specific fine-tuning. We also perform a scaling analysis, showing that by relaxing the constraints on the classical simulated annealing algorithms considered in Ebadi et al., our implementation is competitive with the quantum algorithms. Conversely, instances with larger connectivity or less structure are shown to display a time-to-solution potentially orders of magnitudes larger. Based on these results, we propose protocols to systematically tune problem hardness, motivating experiments with Rydberg atom arrays on instances orders of magnitude harder (for established classical solvers) than previously studied.

Figure 1

Schematic illustration of the problem. (a) We consider unit-disk graphs with nodes arranged on a two-dimensional square lattice with lattice spacing $a$ and filling fraction $\varrho \sim 80\%$, and edges connecting all pairs of nodes within a unit distance (illustrated by the circle). For $\sqrt{2}a\le R_b<2a$ (as considered here), nodes are connected to nearest and next-nearest neighbors resulting in a (quasiplanar) Union-Jack pattern with maximum degree $d_{\max }=8$. (b) Our goal is to solve the MIS problem on this family of instances (as depicted here with nodes colored in red in the right panel) and assess the hardness thereof using both exact and heuristic algorithms.

Figure 2

Schematic illustration of the (exact) sweeping line algorithm (SLA) as applied to the MIS-UD problem. (a) SLA proceeds by sweeping a fictitious line across the graph and tracking all potentially optimal MIS configurations on this boundary, efficiently exploiting the quasiplanar structure of the UD graph. Processed nodes are shown in gray, boundary nodes in blue, and unprocessed nodes in yellow. (b) The light blue node from (a) is added to the boundary, while the bottom left blue node is dropped (as a result of not having any more connections to unprocessed nodes).

Figure 3

Schematic illustration of the heuristic simulated annealing (SA) solver. The original configuration (a) is overlaid with connectivity statistics in (b): nodes in the set (red), nodes without marked neighbors (white), and nodes with a single neighbor (blue). Potential moves are (i) removal of (red) nodes currently in the set, (ii) addition of currently white nodes to the set, and (iii) swapping a blue node with its red neighbor. Gray nodes have more than one adjacent node in the set and no valid moves. From (b) to (c), one node is added to the set and the statistics are updated accordingly. From (c) to (d), one blue node is swapped with its adjacent red node.

Figure 4

Time-to-solution (TTS) for the exact solvers. (a) TTS for the exact SLA solver as a function of system size $N$. For every system size $N$, 1000 random UD instances with $\varrho =0.8$ have been considered. The data fit reasonably to $\mathrm{TTS}\left(N\right)\approx cN{\varphi }^{\sqrt{N}}$, where the basis of $\varphi \approx 1.62$ is the theoretical expectation for the Fibonacci sequence. At larger system sizes ($N>500$), high memory usage causes slower access times (cache misses), resulting in a substantially larger prefactor $c^{\prime }$. (b) TTS for the B&B solver as a function of system size $N$. For every system size $N$, 1000 random UD instances have been considered; see Appendix pp1-s3 for a box plot description. Problems with hundreds (thousands) of nodes can be solved to optimality in subsecond (minute) timescales. The solid line is the linear regression over instances whose TTS are in the top highest $2\%$. The linear regression minimizes the residual sum of squares of log(TTS).

Figure 5

Time required to reach 99% success probability (${\mathrm{TTS}}_{99}$) for the heuristic SA solver as a function of system size $N$ (i.e., how long the solver should run for a 99% chance of finding the optimal solution). For every system size $N$, 1000 random UD instances at $\varrho =0.8$ filling have been considered; see Appendix pp1-s3 for a box plot description. The solid line is the linear regression over instances whose TTS are in the top highest $2\%$. The linear regression minimizes the residual sum of squares of log(TTS).

Figure 6

(a) Time-to-solution (TTS) for the exact SLA solver as a function of the hardness parameter $\mathbb{H}$. Virtually no dependence on $\mathbb{H}$ is observed, showing that TTS is fully determined by the system size $N\sim L^2$. (b) Conversely, for the Markov-chain based SA solver, ${\mathrm{TTS}}_{99}$ shows a strong correlation with the hardness parameter $\mathbb{H}$, as expected.

Figure 7

Estimated success probability $P_{\mathrm{MIS}}$ for the heuristic SA solver as a function of the hardness parameter $\mathbb{H}$. Here we plot $-log(1-P_{\mathrm{MIS}})$ for SA with a fixed depth of 32, for UD graphs selected from the top two percentile of hardness parameter $\mathbb{H}$ for each system size $L=13,\cdots ,33$. Hollow points represent our SA implementation with bias moves in valid configuration space. Power-law fits to the form $\sim {\mathbb{H}}^{-\alpha }$ are used to extract scaling performance with graph hardness $\mathbb{H}$.

Figure 8

(a) Hardness transition as a function of the disk radius (in units of the lattice spacing) $r=R_b/a$, as given by the time-to-solution (TTS) for the B&B solver, shown here for system size $L=21$ and density $\varrho =0.8$ (i.e., $N\approx 350$), with 100 random seeds per radius. (b) TTS as a function of system size $N=\varrho L^2$ for $r=\sqrt{2}$ (blue) and $r=3$ (green), the latter referring to the pronounced peak observed in (a). The solid lines are a linear regression fit over 100 instances with TTS in the highest $10\%$, with corresponding $R^2$ values of 0.87 and 0.98 for $r=\sqrt{2}$ and $r=3$, respectively. Instances with $r=3$ appear to be much harder than those with $r=\sqrt{2}$. See Appendix pp1-s3 for a box plot description.

Figure 9

Hardness transition from unit-disk (UD) to random Erdős-Rényi (ER) graphs. Time-to-solution (TTS) for MIS as a function of the fraction $\varepsilon$ of edges rewired, with 150 random seeds. Starting from Union-Jack-type UD graphs (left), edges are randomly selected and rewired, thereby gradually breaking the UD connectivity, and ultimately generating random ER graphs (right). While the original UD graphs can be solved to optimality in $\sim {10}^{-2}\mathrm{s}$, comparable ER graphs (with the same number of nodes and edges) display a TTS orders of magnitude larger. The red line and the two shaded areas refer to the median TTS over 500 instances for MIS on random ER graphs, the TTS among the $25\%$ and $75\%$ and the minimum and maximum whiskers, respectively. Numerical parameters: $L=21$ and $\varrho =80\%$.

Figure 10

Coefficient $\alpha$ extracted from time-to-solution (TTS) scaling with system size as $\mathrm{TTS}=O\left(2^{\alpha N}\right)$ for the B&B solver as a function of the lattice filling $\varrho$. In the main text, we have focused on $\varrho =0.8$. All results have been taken from the top 2% TTS instances over 1000 random instances. Error bars correspond to the $95\%$ confidence interval.

Figure 11

Correlation analysis for B&B solver. Results are shown for the hardest $2\%$ according to the hardness parameter $\mathbb{H}$ on 1000 instances for each (odd) problem size from $L=7$ to 35. (a) Scatter plot of hardness $\mathbb{H}$ and TTS for B&B solver. The corresponding Pearson correlation amounts to 0.81. (b) Scatter plot of the residuals from the linear regression of $\mathbb{H}$ and TTS with system size $N$. Here the Pearson correlation is found to be 0.01.

Figure 12

Time-to-solution (TTS) for the SLA solver as a function of the disk radius $r$ for random UD instances, for system size $L=21$ and density $\varrho =0.8$ (i.e., $N\approx 350$), with 100 random seeds per radius. Similar to our results for our B&B solver, we observe distinct peaks at $r=2,3,4$.

Figure 13

Hardness parameter $\mathbb{H}$ as a function of the radius of the unit-disk graph (UDG), for system size $L=21$ and density $\varrho =0.8$ (i.e., $N\approx 350$), with 1000 random seeds per radius.