Exponential Speed-ups for Structured Goemans-Williamson relaxations via Quantum Gibbs States and Pauli Sparsity

TL;DR

The paper identifies a class of QUBO instances for which the Goemans–Williamson SDP can be solved exponentially faster using matrix multiplicative weight methods together with quantum Gibbs-state simulations, provided the cost matrix is Pauli-sparse with a small diagonal group. By relaxing constraints to a carefully chosen subset and exploiting locality through Pauli decompositions, the authors show polylogarithmic-time solvability in the problem dimension $D=2^n$ under certain structural conditions, and they develop practical rounding schemes to extract high-quality QUBO energies via Monte Carlo methods. They demonstrate end-to-end performance on very large instances (up to $D=2^{50}$) using tensor networks and Kronecker graphs as testbeds, yielding energies within about $0.15\\%$ of the GW optimum on challenging scales. The work bridges convex optimization, quantum many-body physics, and combinatorial algorithms, offering both rigorous exponential speedups for structured SDPs and heuristic methods with broad applicability. It also highlights open questions about real-world applicability, optimal constraint selection, and the potential expansion to broader graph families and higher local dimensions.

Abstract

Quadratic Unconstrained Binary Optimization (QUBO) problems are prevalent in various applications and are known to be NP-hard. The seminal work of Goemans and Williamson introduced a semidefinite programming (SDP) relaxation for such problems, solvable in polynomial time that upper bounds the optimal value. Their approach also enables randomized rounding techniques to obtain feasible solutions with provable performance guarantees. In this work, we identify instances of QUBO problems where matrix multiplicative weight methods lead to quantum and quantum-inspired algorithms that approximate the Goemans-Williamson SDP exponentially faster than existing methods, achieving polylogarithmic time complexity relative to the problem dimension. This speedup is attainable under the assumption that the QUBO cost matrix is sparse when expressed as a linear combination of Pauli strings satisfying certain algebraic constraints, and leverages efficient quantum and classical simulation results for quantum Gibbs states. We demonstrate how to verify these conditions efficiently given the decomposition. Additionally, we explore heuristic methods for randomized rounding procedures and extract the energy of a feasible point of the QUBO in polylogarithmic time. While the practical relevance of instances where our methods excel remains to be fully established, we propose heuristic algorithms with broader applicability and identify Kronecker graphs as a promising class for applying our techniques. We conduct numerical experiments to benchmark our methods. Notably, by utilizing tensor network methods, we solve an SDP with $D = 2^{50}$ variables and extract a feasible point which is certifiably within $0.15\%$ of the optimum of the QUBO through our approach on a desktop, reaching dimensions millions of times larger than those handled by existing SDP or QUBO solvers, whether heuristic or rigorous.

Figures (6)

• Figure 1: Performance of different solvers on $C$ given by a non-trivial 1D commuting Hamiltonian [see Eq. \ref{['eq:1d_hamiltonian']}]: (a) GW SDP and rounded solutions (see \ref{['app:details_numerics_tensor']} for technique details of calculating the SDP and performing the randomized rounding with tensor network approach). (b) Ratio between the rounded and SDP solutions from HU solver.$\operatorname{GW}(C,\dot{\mathcal{D}}_C,\epsilon)$ denotes the SDP solution obtained from the HU solver. RR $\operatorname{GW}(C, \dot{\mathcal{D}}_C,\epsilon)$ refers to the mean objective value over rounded strings sampled from the $\operatorname{GW}(C, \dot{\mathcal{D}}_C,\epsilon)$, corresponding to a feasible, approximate solution of the associated $\operatorname{QUBO}$ problem. The objective value $\operatorname{QUBO}(C)$, see \ref{['eq:GW_relaxation_normalized']}] of the exact solution is unknown, but it is guaranteed to lie within RR $\operatorname{GW}(C, \dot{\mathcal{D}}_C,\epsilon)$ and $\operatorname{GW}(C, \dot{\mathcal{D}}_C,\epsilon)$, up to additive error $\epsilon$. The Monte Carlo sampling error for the randomized rounding is $\epsilon=0.5$---see \ref{['equ:complexity_bound_density']} see for details. In turn, $\operatorname{GW}(C,\varnothing ,\epsilon)$ represents the solution to the SDP spectral relaxation. SketchyCGAL and CVXPY denote the $\operatorname{GW}(C)$ SDP solutions obtained via SketchyCGAL Yurtsever2021 and the standard CVXPY SDP Python library diamond2016cvxpy. The maximum problem sizes they can handle are $D=2^{14}$ and $D=2^9$, respectively, due to memory limitations. Additionally, we plot the objective values $\operatorname{QUBO}(C)$ obtained using the standard Tabu search Python library (dwave-tabu) with 50 samples stably up to $D=2^{10}$. In the numerics, we observe that Tabu starts crashing frequently starting at $D=2^{11}$. Results are averaged over 100 randomly generated problem instances, as we observe a high degree of concentration for the value of the random instances. The error bars, except for those in $\operatorname{GW}(C,\dot{\mathcal{D}}_C,\epsilon)$, represent the standard error of the mean across all instances. The error bar in $\operatorname{GW}(C,\dot{\mathcal{D}}_C,\epsilon)$ denotes the error tolerance introduced in HU, which dominates the numerical errors. All numerical experiments run on a single, standard desktop computer with specifications in \ref{['app:details_numerics']}.
• Figure 2: Lower bounds to the tree- and rank-widths of the instances considered for Fig. \ref{['fig:sdp_1d_cluster']} computed through convex relaxations Cygan2015 Both lower bounds feature of a super-linear growth in the number of vertices (dimension $D$ of $C$). In addition, the genus of the graphs (not shown) is observed to grow even significantly faster. This implies that the known algorithms able to exploit structures associated to these complexity parameters cannot solve the corresponding QUBO instances in time less than super-polynomial in $D$ (see main text). Hence, this sanity check is consistent with the instances defined by Eq. \ref{['eq:1d_hamiltonian']} being computationally non-trivial. Note that these parameters only depend on the graph on which these instances are defined on, not on the particular choice of the weights $a_i$.
• Figure 3: Performance of different solvers on $\tilde{C}$ sampled from a Kronecker graph $A^{\otimes k}$ for $k=1,2,3,4,5,6$ [\ref{['eq:kronecker_graph']}]: (a) GW SDP and rounded solutions(see \ref{['app:details_numerics_sparse']} for technique details of calculating the SDP and performing the randomized rounding(label with prefix RR) with sparse matrix-vector approach). (b) Ratio between the rounded and SDP solutions from HU solver. (c) Number of constraints used in the HU solver. The Pauli-sparse approximation $\tilde{C}$ is generated using the sampling method in \ref{['thm:pauli_sparsification']} with the number of samples determined by \ref{['eq:exact_sample_kronecker']}, and is then renormalized as $\tilde{C}/\|\tilde{C}\|_{P,l_1}$. The trace and the rounding are estimated using $50$ random vectors with the sparse matrix-vector multiplication method described in \ref{['subsec:poly_speedups']}. $\operatorname{GW}(C,\mathcal{D}_{\tilde{C},2},\epsilon)$ and $\operatorname{GW}(C,\mathcal{D}_{\tilde{C},3},\epsilon)$ denote the SDP solutions obtained by the HU solver with constraint sets derived from the Krylov heuristics method (\ref{['sec:approximate_constraint_set']}) up to second and third order, respectively, and the corresponding error bars in (a) indicate the maximum error introduced within the HU framework. RR $\operatorname{GW}(C,\mathcal{D}_{\tilde{C},2},\epsilon)$ and RR $\operatorname{GW}(C,\mathcal{D}_{\tilde{C},3},\epsilon)$ refer to the mean objective value over rounded strings sampled from $\operatorname{GW}(C,\mathcal{D}_{\tilde{C},2},\epsilon)$ and $\operatorname{GW}(C,\mathcal{D}_{\tilde{C},3},\epsilon)$, respectively, corresponding to feasible approximate solutions of the associated $\operatorname{QUBO}$ problem. SketchyCGAL and CVXPY denote the $\operatorname{GW}(C)$ SDP solutions obtained via SketchyCGAL Yurtsever2021 and the standard CVXPY SDP Python library diamond2016cvxpy. The maximum problem sizes they can handle are $D=2^{12}$ and $D=2^8$, respectively, due to memory limitations. Additionally, we plot the objective values $\operatorname{QUBO}(C)$ obtained using the standard Tabu search Python library (dwave-tabu) with 300 samples. The results show a clear drop in performance at $D=2^{12}$, with error bars representing sampling variability in the Tabu procedure. All numerical experiments run on a single, standard desktop computer with specifications in \ref{['app:details_numerics']}.
• Figure 4: Parameterized tensors associated with the partition function of the 1D Hamiltonians $\operatorname{GW}(C,S,\epsilon)$[\ref{['eq:partition_function_main']}], together with the additional tensor nodes used for randomized rounding. The explicit mathematical expressions corresponding to each tensor are displayed beneath the respective nodes. Square nodes denote tensors acting on multiple qubits, triangular nodes represent single-qubit vector tensors, and circular nodes correspond to single-qubit unitaries. Colors and text labels are used to differentiate tensors beyond shape alone. The index $I^\mathrm{i}_{q_i}$ indicates the input leg associated with the $q_i$-th qubit, while $I^\mathrm{o}_{q_i}$ denotes the corresponding output leg.
• Figure 5: Tensor network representation of $\exp\!\left(H(\lambda,C,S)\right)$. The network is built by sequentially contracting the tensors associated with the terms in $\exp(-\lambda_{C} \frac{C}{\|C\|})$ along their shared indices (large pink square), followed by contraction with the tensors corresponding to the constraint Hamiltonians (large grey square). External legs $I_1^{\mathrm{i}}, I_2^{\mathrm{i}}, \dots, I_n^{\mathrm{i}}$ and $I_1^{\mathrm{o}}, I_2^{\mathrm{o}}, \dots, I_n^{\mathrm{o}}$ denote input and output indices for each qubit. The trace of the operator is obtained by contracting each input index $I_i^{\mathrm{i}}$ with the output index $I_i^{\mathrm{o}}$ iteratively.

Theorems & Definitions (65)

• Definition 1: Diagonal group of $C$
• Definition 2: Pauli-sparse
• Remark 1
• Remark 2
• Remark 3
• Remark 4
• Definition 3: Pauli $\ell_1$ norm
• Definition 4: Relaxed Goemans-Williamson
• Theorem 5: Reducing the number of constraints, informal
• Corollary 6: Condition for spectral relaxation, informal