Solving quadratic binary optimization problems using quantum SDP methods: Non-asymptotic running time analysis

TL;DR

An analysis of the non-asymptotic resource requirements of the quantum SDP solvers proposed, which includes a version with adaptive step-sizes, an improved detection criterion for infeasible instances, and a more efficient rounding procedure.

Abstract

Quantum computers can solve semidefinite programs (SDPs) using resources that scale better than state-of-the-art classical methods as a function of the problem dimension. At the same time, the known quantum algorithms scale very unfavorably in the precision, which makes it non-trivial to find applications for which the quantum methods are well-suited. Arguably, precision is less crucial for SDP relaxations of combinatorial optimization problems (such as the Goemans-Williamson algorithm), because these include a final rounding step that maps SDP solutions to binary variables. With this in mind, Brandão, França, and Kueng have proposed to use quantum SDP solvers in order to achieve an end-to-end speed-up for obtaining approximate solutions to combinatorial optimization problems. They did indeed succeed in identifying an algorithm that realizes a polynomial quantum advantage in terms of its asymptotic running time. However, asymptotic results say little about the problem sizes for which advantages manifest. Here, we present an analysis of the non-asymptotic resource requirements of this algorithm. The work consists of two parts. First, we optimize the original algorithm with a particular emphasis on performance for realistic problem instances. In particular, we formulate a version with adaptive step-sizes, an improved detection criterion for infeasible instances, and a more efficient rounding procedure. In a second step, we benchmark both the classical and the quantum version of the algorithm. The benchmarks did not identify a regime where even the optimized quantum algorithm would beat standard classical approaches for input sizes that can be realistically solved at all. In the absence of further significant improvements, these algorithms therefore fall into a category sometimes called galactic: Unbeaten in their asymptotic scaling behavior, but not practical for realistic problems.

Figures (7)

• Figure 1: Illustration of the Hamiltonian Updates algorithm. The circle depicts the space of trace-one psd matrices, with $\rho_0$ lying in the center. The feasible region is shown in dark red, the $\epsilon$-feasible region is marked light red. Each graphic shows a single iteration of HU: At the start of each update, a matrix $P$ is calculated that defines a hyperplane (shown in green) separating the current $\rho$ and the feasible region. By updating $\rho_H\rightarrow\rho_{H+\lambda P}$ (i.e. penalizing infeasible directions), $\rho$ moves towards the hyperplane. For simplicity this is depicted by a straight line. This is generally not the case, as $\rho$ depends non-linearly on $H$. The procedure ends when $\rho$ enters the $\epsilon$-feasible region (i.e. all constraints are fulfilled up to precision $\epsilon$ as defined in \ref{['eq:approx_feasibility_constraint']}).
• Figure 2: Behavior of the improved HU algorithm on an instance used with parameters $n=1024$, $s=16$ and $\epsilon=0.001$, generated as described in Sec. \ref{['sec:numerics_improvements']}. The algorithm terminated after finding an $\epsilon$-feasible solution for $\gamma:=\gamma^\star+\mathrm{offset}$, or once the free energy became positive. Here, $\gamma^\star$ is the optimal objective value determined by an SDP solver ODonoghue2016. The error bars show the upper and lower quartiles. Upper panel: Number of iterations required until convergence as a function of the $\gamma$-offset. The improved termination criterion comes into play on the right tail of the curve. We observe that the number of iterations required to certify infeasibility goes down the further the problem specification is from a feasible one. (This contrasts to the fixed termination criterion used in Brandao2022fasterquantum). Lower panel: The lower bound on the free energy at time of termination. We observe numerically that the bound increases as $\gamma$ gets closer to the optimal value. Whether this effect can be exploited for algorithmic improvements is a question we leave open.
• Figure 3: Number of iterations for varying values of the momentum hyperparameter $\beta$. The two different approaches for the diagonal update are compared: The original $\ell_1$ norm based $P_d^{\ell_1}$ (blue) and the new $\ell_2$ norm based $P_d^{\ell_2}$ (orange). The instances have dimension $n=1024$, sparsity $s=16$ and precision $\epsilon=0.001$, and use the optimal SDP solution $\gamma^\star$ as the target objective value. The error bars show the upper and lower quartiles.
• Figure 4: Total number of iterations over a complete binary search as a function of the precision $\epsilon$. The fit is given by $f(\epsilon)=6.48\times 10^{-3} \epsilon^{-2.02}$ with a 95% confidence interval of $(-2.38, -1.73)$ for the exponent of $\epsilon$. Therefore, the experimentally observed scaling exponent is similar to the theoretical bound of $-2$ (however, the experimentally observed prefactor is significantly better). The error bars show the upper and lower quartiles. The gray area shows the 95% confidence interval of the fit.
• Figure 5: Behavior of the objective values $x^TC x / n$ after rounding. Results are obtained by performing a binary search over HU instances for different precision parameters $\epsilon$ and applying randomized rounding. The results are averaged over 20 cost matrices with dimension $n=1024$ and sparsity $s=16$. The error bars show the upper and lower quartiles. Left panel: Comparison of average objective value (blue) and the maximum value seen in 1000 roundings (orange). Right panel: Difference $\nu=(x_\mathrm{opt}^TC x_\mathrm{opt}- x_\epsilon^TC x_\epsilon)/n$, where $x_\mathrm{opt}$ and $x_\epsilon$ have been obtained, respectively, by applying the rounding procedure to an optimal SDP solution, and an $\epsilon$-precise HU solution. The fit for the average value is given by $f_{\nu}(\epsilon)=0.32 \epsilon^{0.79}$ with a 95% confidence interval of $(0.71,0.88)$ for the exponent. The colored areas show the 95% confidence interval of the fits.

Theorems & Definitions (30)

• Theorem 1: Brandao2022fasterquantum
• Lemma 1
• Theorem 2
• Theorem 3
• Theorem 4
• Lemma 5: vanApeldoorn2020quantumsdpsolvers
• Lemma 5
• proof
• Lemma 6
• proof