Simulated Annealing for Quadratic and Higher-Order Unconstrained Integer Optimization

TL;DR

The paper tackles solving polynomial unconstrained optimization with integer variables (QUIO/HUIO) by applying simulated annealing directly, avoiding the combinatorial blow-up from binary encoding to QUBO/HUBO. It introduces an efficient SA framework that stores and updates energy-difference coefficients, enabling fast $O(1)$ updates for QUIO and $O(|\Omega_k|)$ updates for HUIO, and adds an optimal-transition Metropolis mechanism to cope with wide variable ranges. Through comprehensive benchmarks on fully-connected, multilinear, and random-interaction models, the authors demonstrate that direct QUIO/HUIO SA typically achieves better solution quality than QUBO-based SA, with the optimal-transition variant providing notable advantages for short annealing times. The approach is implemented in the OpenJij library, offering a practical tool for scalable integer polynomial optimization in real-world applications.

Abstract

Simulated annealing (SA) is a key algorithm for solving combinatorial optimization problems, which model numerous real-world systems. While SA is commonly used to solve quadratic unconstrained binary optimization (QUBO) problems, many practical problems are more naturally formulated using integer variables. We therefore study quadratic and higher-order unconstrained integer optimization (QUIO and HUIO) problems, which generalize QUBO by allowing integer-valued variables and higher-order interactions. Conventional approaches often convert these problems into QUBO formulations through binary encoding and the reduction of higher-order terms. Such conversions, however, greatly increase the number of variables and interactions, resulting in long computation times and, for large-scale problems, even making the conversion itself a dominant computational bottleneck. To overcome this limitation, we propose an efficient framework that directly applies SA to QUIO and HUIO problems without converting them into QUBO. Within this framework, we introduce an optimal-transition Metropolis method, designed to improve efficiency when the variable range is wide, and evaluate its performance alongside the conventional Metropolis and heat bath methods. Numerical experiments demonstrate that the proposed direct approach achieves higher efficiency and better solution quality than the conventional QUBO-based formulation and reveal the practical advantages of the optimal-transition Metropolis method. The algorithm developed in this study is available as part of the open-source library OpenJij, which provides a Python front end with a C++ backend.

Figures (8)

• Figure 1: Annealing time as a function of variable range $2u + 1$ for (a) $E^{\text{FC}}_{p=1,u}(\bm{z})$ and (b) $E^{\text{FC}}_{p=2,u}(\bm{z})$. Each data point represents the mean annealing time over 100 samples, with 1000 sweeps and $N = 100$ variables. Error bars indicate the standard error.
• Figure 2: Annealing time as a function of $p$ for (a) $E^{\text{FC}}_{p,u=1}(\bm{z})$ and (b) $E^{\text{ML}}_{p,u=10}(\bm{z})$. Each data point represents the mean annealing time over 100 samples, with 1000 sweeps and $N = 100$ variables. Error bars indicate the standard error but are too small to be visible.
• Figure 3: Temperature $T$ dependence of the mean energy for (a) $E^{\text{FC}}_{p=2,u=1}(\bm{z})$, (b) $E^{\text{FC}}_{p=2,u=100}(\bm{z})$, (c) $E^{\text{RI}}_{p=4,u=1}(\bm{z})$, and (d) $E^{\text{RI}}_{p=4,u=100}(\bm{z})$. Each data point represents the mean energy at temperature $T$ during SA, averaged over 1000 samples, and the number of sweeps is 1000. The number of variables is $N = 100$, and for $E^{\text{RI}}_{p,u}(\bm{z})$, the number of interactions is 1000. The black dotted line for (a) and (b) indicates the ground-state energy, $-5050$.
• Figure 4: Temperature $T$ dependence of the mean energy for $E^{\text{ML}}_{p=30,u=1}(\bm{z})$. Each data point represents the mean energy at temperature $T$ during SA, averaged over 1000 samples, and the number of sweeps is (a) 1000 and (b) 100000. The number of variables is $N = 100$, and the black dotted line indicates the ground-state energy, $-100$.
• Figure 5: Sweep dependence of the mean energy for the random interaction model [Eq. (\ref{['model-RI']})] with (a) $p=2,u=1$ and (b) $p=2,u=100$. Each data point represents the mean energy averaged over 1000 samples. Error bars indicate the standard error but are too small to be visible. The number of variables is $N = 100$ and the number of interactions is 1000.