Simulated bifurcation for higher-order cost functions

Special purpose hardware solvers for the Ising problem Barahona1982 have been applied to combinatorial optimization problems Lucas2014 . Such hardware solvers are referred to as Ising machines Mohseni2022a . Ising machines with analogue computers have attracted attention, such as quantum annealers Kadowaki1998 with superconducting circuits Johnson2011 and coherent Ising machines with pulse lasers Marandi2014 . Also, Ising machines with digital processors have been developed based on, e.g., simulated annealing (SA) Kirkpatrick1983 ; Yamaoka2016 ; Aramon2019 , simulated bifurcation (SB) Goto2019 , and other algorithms Leleu2021a .

SB is a heuristic algorithm based on numerical simulation of a classical counterpart of a quantum bifurcation machine, which is a quantum annealer using nonlinear oscillators Goto2016 ; Puri2017 ; Goto2019a ; Kanao2021 . SB allows us to simultaneously update variables, and thus can easily be accelerated by parallel processing with, e.g., graphics processing units (GPUs) and field-programmable gate arrays (FPGAs) Goto2019 ; Tatsumura2019 ; Zou2020 ; Tatsumura2021 . Two improved variants of SB have been reported: ballistic and discrete SBs (bSB and dSB) Goto2021 . In addition, bSB and dSB with thermal fluctuations have been studied Kanao2022a . SB has been benchmarked against other Ising machines Oshiyama2022 ; Kowalsky2022 . Besides, SB has been applied to financial portfolio optimizations Steinhauer2020a ; Cohen2020 , financial trading machines Tatsumura2020 , traveling salesman problems Zhang2021c , clustering with a hybrid method Matsumoto2022 , and wireless communications Zhang2022c .

In the Ising problem, the aim is to minimize cost functions written in second-order polynomials of spin variables, for which Ising machines have been designed. However, some combinatorial optimization problems are naturally expressed by cost functions with higher-order polynomials than second order (higher-order cost functions). To apply Ising machines to such problems, higher-order cost functions have been transformed to second-order ones Kolmogorov2004 ; Ishikawa2011 ; Xia2018 ; Dattani2019 ; Fujisaki2022a , though such transformations usually require additional spin variables and thus increase computational costs.

In this work, we extend SB to apply it directly to a problem with higher-order cost functions. Since SB is based on simulation, this approach can easily be implemented. To benchmark the higher-order SB, we use a three-regular three-XORSAT (3R3X) problem Hen2019 ; Bernaschi2021 ; Bellitti2021 ; Kowalsky2022 , because this problem yields third-order cost functions and its solution can be known in advance. We compare the third-order SB with that transformed to second order. As a result, we find that third-order bSB and dSB (3bSB and 3dSB) show higher performance than the second-order ones (2bSB and 2dSB). Surprisingly, 3bSB performs better than 3dSB, while 2dSB does better than 2bSB as in the previous study Goto2021 . We further find that 3bSB can be superior to SA applied directly to the third-order cost functions. These results suggest that the higher-order SB will offer another useful approach to special purpose hardware solvers for combinatorial optimization problems.

A third-order cost function with $N$ spin variables $s_{i}=\pm 1$ can generally be written as

$\displaystyle E(\bm{s})$

$\displaystyle=$

$\displaystyle-\!\sum_{i}\!h_{i}s_{i}-\!\sum_{i,j}\!J_{ij}s_{i}s_{j}-\!\sum_{i,j,k}\!K_{ijk}s_{i}s_{j}s_{k},$

(1)

where $\bm{s}$ is a vector consisting of $s_{i}$, and $h_{i},J_{ij}$, and $K_{ijk}$ are coefficients. Since $s_{i}^{2}=1$, we assume without loss of generality that $J_{ii}=0$, and $K_{ijk}=0$ if two of $i,j$, and $k$ are equal.

We here extend second-order bSB and dSB to third-order ones. The third-order bSB and dSB are based on the following equations of motion for positions $x_{i}$ and momenta $y_{i}$ corresponding to $s_{i}$:

	$\displaystyle\dot{y}_{i}$	$\displaystyle=$	$\displaystyle-\left[a_{0}-a(t)\right]x_{i}+cf_{i},$		(2)
	$\displaystyle\dot{x}_{i}$	$\displaystyle=$	$\displaystyle a_{0}y_{i},$		(3)

where the dots denote time derivatives, $a_{0}$ is a constant parameter, $a(t)$ is a bifurcation parameter increased from 0 to $a_{0}$ with the time $t$, and $c$ is a normalization factor determined from forces $f_{i}$. The forces $f_{i}$ are given by

	$\displaystyle f_{i}$	$\displaystyle=$	$\displaystyle-\frac{\partial E(\bm{x})}{\partial x_{i}}=G_{i}(\bm{x})\quad\text{(bSB)},$		(4)
	$\displaystyle f_{i}$	$\displaystyle=$	$\displaystyle-\left.\frac{\partial E(\bm{x})}{\partial x_{i}}\right|_{\bm{x}=\bm{s}}=G_{i}(\bm{s})\quad\text{(dSB)}.$		(5)

For 3bSB, $f_{i}$ are the gradients $G_{i}(\bm{x})$ of $E(\bm{x})$, while for 3dSB, $\bm{x}$ in the gradients is substituted by $\bm{s}$, which is a vector of the signs of $x_{i}$, $s_{i}={\rm sgn}\left(x_{i}\right)$. Both 3bSB and 3dSB assume perfectly inelastic walls at $x_{i}=\pm 1$ Goto2021 , that is, if $\left|x_{i}\right|>1$, $x_{i}$ and $y_{i}$ are set to ${\rm sgn}\left(x_{i}\right)$ and $0$, respectively. Equations (2) and (3) are solved with the symplectic Euler method Leimkuhler2004 , which discretizes the time with an interval $\Delta t$ and computes time evolutions of $x_{i}$ and $y_{i}$ step by step. The signs of $x_{i}$ at the final time give a solution. The solution at least corresponds to a local minimum of $E(\bm{s})$ as for 2bSB and 2dSB Goto2021 , because Eq. (2) becomes $\dot{y}_{i}=c_{0}f_{i}$ at the final time. Note that when all $K_{ijk}$ are zero in Eq. (1), Eqs. (2)-(5) reproduce 2bSB and 2dSB in Ref. Goto2021, . We thus focus on the third term in Eq. (1) in the following.

To reduce the computational cost of calculating $f_{i}$, we utilize sparse connectivity of the 3R3X problem, that is, the number, $M$, of nonzero $K_{ijk}$ is much smaller than the total number of $K_{ijk}$. We expect that sparse connectivities often appear in applications and therefore the following method will be useful. By denoting nonzero $K_{ijk}$ by $K_{m}$ $(m=1,2,\cdots,M)$, $E(\bm{s})$ is expressed as

$\displaystyle E(\bm{s})$

$\displaystyle=$

$\displaystyle-\sum_{m=1}^{M}K_{m}s_{v_{m1}}s_{v_{m2}}s_{v_{m3}},$

(6)

where $v_{mn}$ with $n=1,2,3$ represent the indices of spin variables included in the $m$th term. To further reduce the time required for evaluating the third-order SB, we calculate the gradients by

	$\displaystyle G_{i}(\bm{x})$	$\displaystyle=$	$\displaystyle\frac{1}{x_{i}+\epsilon}\sum_{m=1}^{M}a_{im}K_{m}x_{v_{m1}}x_{v_{m2}}x_{v_{m3}}\quad\text{(bSB)},$		(7)
	$\displaystyle G_{i}(\bm{s})$	$\displaystyle=$	$\displaystyle s_{i}\sum_{m=1}^{M}a_{im}K_{m}s_{v_{m1}}s_{v_{m2}}s_{v_{m3}}\quad\text{(dSB)},$		(8)
	$\displaystyle a_{im}$	$\displaystyle=$	$\displaystyle\delta_{i,v_{m1}}+\delta_{i,v_{m2}}+\delta_{i,v_{m3}},$		(9)

where $\delta_{i,j}$ is 1 if $i=j$ and otherwise 0. Equation (8) exactly gives the derivative in Eq. (5), because $a_{im}$ pick up terms including $s_{i}$ from Eq. (6), and these $s_{i}$ are canceled by multiplying $s_{i}$ in Eq. (8) with $s_{i}^{2}=1$. This method is valid since Eq. (6) is linear with respect to each $s_{i}$. Similarly, Eq. (7) yields the same value as Eq. (4) if $x_{i}\neq 0$ and $\epsilon=0$. To stabilize calculations, we set $\epsilon=10^{-14}$. We empirically found that Eqs. (7) and (8) could reduce the evaluation time compared with direct differentiation, which may be because these equations gather three terms derived from $x_{v_{m1}}x_{v_{m2}}x_{v_{m3}}$ together and thus reduce memory accesses for $x_{v_{mn}}$ and $s_{v_{mn}}$.

The 3R3X problem is mapped to a minimization problem of Eq. (6) with the following $v_{mn}$ and $K_{m}$ Hen2019 . Here, $M$ equals $N$. For each of $n=1,2,3$, the vector of indices $v_{mn}$ is a permutation of $1,2,\cdots,N$. For each $m$, three indices $v_{m1}$, $v_{m2}$, and $v_{m3}$ are different from each other. For different $m$, the combinations of $v_{m1}$, $v_{m2}$, and $v_{m3}$ are different. $K_{m}$ is given by $(-1)^{b_{m}}$ with $b_{m}=0,1$. Thus $E(\bm{s})\geq-N$ holds.

For a 3R3X instance, a solution $\bm{s}_{0}$ minimizing $E(\bm{s})$ can be obtained with linear algebra Hen2019 , that is, $\bm{s}_{0}$ is related to a solution, $\xi_{i}=0,1$, of linear equations,

$\displaystyle\sum_{i=1}^{N}a_{mi}\xi_{i}$

$\displaystyle=$

$\displaystyle b_{m},\quad\mathrm{mod}\quad 2,$

(10)

by $s_{i0}=(-1)^{\xi_{i}}$. However, depending on $a_{mi}$ and $b_{m}$, Eq. (10) may have no solution. To ensure that Eq. (10) has at least one solution, we instead give $a_{mi}$ and $\xi_{i}$, and determine $b_{m}$ with Eq. (10). Then an exact solution $\bm{s}_{0}$ is indicated by $E(\bm{s}_{0})=-N$.

A third-order cost function of the 3R3X problem can be decomposed to second-order one by mapping a term $-(-1)^{b}s_{1}s_{2}s_{3}$ to

	$\displaystyle h\left(s_{1}+s_{2}+s_{3}\right)+\tilde{h}\tilde{s}$
	$\displaystyle+J\left(s_{1}s_{2}+s_{2}s_{3}+s_{3}s_{1}\right)+\tilde{J}\left(s_{1}+s_{2}+s_{3}\right)\tilde{s},$		(11)

where $\tilde{s}=\pm 1$ is an ancillary spin variable, and $J$, $\tilde{J}$, $h$, and $\tilde{h}$ are coefficients determined so that $s_{1}$, $s_{2}$, and $s_{3}$ minimizing Eq. (11) coincide with ones minimizing $-(-1)^{b}s_{1}s_{2}s_{3}$ Kolmogorov2004 ; Leib2016 ; Hen2019 . We use $J=1/4$, $\tilde{J}=(-1)^{b}/2$, $h=-(-1)^{b}/4$, and $\tilde{h}=-1/2$, which normalize the minimum of Eq. (11) to $-1$. The ancillary spin variables increase the total number of spin variables to $2N$.

We also directly apply SA to third-order cost functions (3SA). In an update in SA Kirkpatrick1983 , one calculates a change in a cost function $\Delta E_{i}$ caused by a flip of $s_{i}$. The flip is accepted if $\beta\Delta E_{i}<-\ln R$, where $\beta$ and $R$ are the inverse temperature increased during annealing and a random number in the interval $(0,1)$, respectively. In a step of SA, we sequentially update spin variables $s_{i}$ Isakov2015 .

In numerical experiments, we measure computational costs of finding exact solutions as follows. We randomly generate 100 3R3X instances for each $N$. For these instances, we run the above mentioned SB and SA many times with different initial conditions or random numbers, and obtain probabilities of finding exact solutions, $P$. Then, using the number of steps in a run, $N_{\mathrm{s}}$, we estimate the number of steps necessary to obtain exact solutions with a probability of 99% (step-to-solution) by Leleu2021a ; Kanao2022a

$\displaystyle S$

$\displaystyle=$

$\displaystyle N_{\mathrm{s}}\frac{\log 0.01}{\log(1-P)}.$

(12)

A smaller $S$ indicates that an algorithm provides an exact solution more efficiently (or faster if computation time for each step is assumed to be the same).

We set the parameters as follows. For SB, $a_{0}=1$, and $a(t)$ is increased linearly. To normalize the term $cf_{i}$ in Eq. (2), we give $c$ by Sakai2020

$\displaystyle c$

$\displaystyle=$

$\displaystyle c_{1}\left(\frac{1}{\nu}\sum_{i=1}^{\nu}f_{i}^{2}\right)^{-1/2},$

(13)

where $c_{1}$ and $\nu$ are a constant parameter and the number of spin variables, respectively ($\nu$ are $2N$ and $N$ for the second- and third-order SBs, respectively). We optimize $N_{\mathrm{s}}$ for each $N$, and $(\Delta t,c_{1})$ for certain large values of $N$, to minimize the median of $S$ for the 100 instances Kowalsky2022 . Initial $x_{i}$ and $y_{i}$ are generated by uniform random numbers in the interval $(-1,1)$. To suppress a statistical error in $P$, the number of runs $N_{\mathrm{r}}$ is set to be large such that $N_{\mathrm{r}}P\gg 1$ Goto2021 , though this condition is not necessarily met for the largest $N$ in this work owing to limitations of computational resources. For SA, $\beta$ is linearly increased from zero to a final value $\beta_{1}$, and $\beta_{1}$ and $N_{\mathrm{s}}$ are optimized as in SB.

We first present typical distributions of $E(\bm{s})$ for a 3R3X instance, for which 3bSB gives $S$ near the median of the 100 instances. Figure 1 shows the distributions by the rates of runs resulting in $E(\bm{s})$. A higher value of the rate means a higher probability of finding solutions with certain accuracy, though here the parameters are optimized to minimize $S$ as mentioned above and the rates can be improved by increasing $N_{\mathrm{s}}$. Interestingly, Fig. 1(a) shows that the rate for small $E(\bm{s})$ is higher for 3bSB than for 3dSB despite smaller $N_{\mathrm{s}}$ for 3bSB. On the other hand, Fig. 1(b) shows that the rate is higher for 2dSB than for 2bSB in accordance with much larger $N_{\mathrm{s}}$ for 2dSB. As shown in Fig. 1(a), the rate for 3SA is between those of 3bSB and 3dSB. These rates and $N_{\mathrm{s}}$ suggest the highest performance of 3bSB.

We next evaluate the results for 100 3R3X instances. Figure 2 shows the medians of $S$ as functions of $N$. All the algorithms exhibit exponential increases of $S$ with $N$. The values of $S$ are smaller for third-order bSB and dSB than for the second-order ones. Among the second-order ones, 2dSB gives smaller $S$ than 2bSB, which is consistent with a reported result that dSB more often reaches exact solutions for the Ising problem Goto2021 . Among the third-order SB, in contrast, 3bSB shows smaller $S$ than 3dSB. The higher performance of 3bSB might be because bSB can rapidly converge to a local minimum of a cost function Goto2021 and a local search algorithm is effective for the 3R3X problem Bernaschi2021 ; Bellitti2021 ; Kowalsky2022 . A comparison between 2bSB and 3bSB suggests that the local search by bSB is especially effective when directly applied to third-order cost functions with fewer spin variables. 3SA results in $S$ between those of 3bSB and 3dSB, indicating that 3bSB is superior to 3SA.

Scaling can be measured by exponents $\alpha$ and $\beta$ such that

$\displaystyle S$

$\displaystyle\simeq$

$\displaystyle 10^{\alpha N+\beta},$

(14)

for large $N$ Kowalsky2022 . A smaller $\alpha$ allows an algorithm to be applied to larger instances. We fit the medians of $S$ by Eq. (14) with the least squares method as the lines in Fig. 2. Here, only the data for intermediate $N$ are used for the fittings, because errors due to finite-size effects and statistical errors are large for small and large $N$, respectively Kowalsky2022 . Table 1 shows $\alpha$ and $\beta$. The value of $\alpha$ for 3bSB is the smallest. $\alpha$ for 3dSB is a little larger than it, and $\alpha$ for 2dSB is significantly larger. $\alpha$ for 2bSB is by far the largest. The value of $\alpha$ for 3SA is larger than those for 3bSB and 3dSB. These $\alpha$ demonstrate that 3bSB scales the best among these algorithms.

Finally, we compare $\alpha$ with the results of Ref. Kowalsky2022, , where various Ising machines are benchmarked using the 3R3X problem expressed by second-order cost functions. Among the Ising machines, Digital Annealer reported the smallest exponent corresponding to $\alpha=0.0370(8)$ Footnote2022 . Compared with this value, $\alpha$ for 3bSB seems to be smaller. $\alpha$ for Simulated Bifurcation Machine (provided by Toshiba Digital Solutions Corporation TDSL ) corresponds to $0.0434(12)$, which coincides with that for 2dSB in this work within the range of errors. This coincidence may support the validity of our evaluation. In Ref. Kowalsky2022, , a quasi-greedy algorithm called SATonGPU Bernaschi2021 ; Bellitti2021 ; Kowalsky2022 is also compared, which is specialized for the 3R3X problem and deals with third-order cost functions directly. As a result, SATonGPU provided the exponent corresponding to $\alpha=0.0341(2)$ Footnote2021 , which is close to an estimated optimal exponent Bernaschi2021 . The value of $\alpha$ for SATonGPU is smaller than $\alpha$ for 3bSB in our work, but the difference is fairly small. This implies that 3bSB can achieve almost optimal scaling even though 3bSB is not specialized for the 3R3X problem.

In summary, we have demonstrated that third-order SB can perform better than second-order one for a 3R3X problem. Although 2dSB provides better performance than 2bSB, 3bSB instead does better than 3dSB. Furthermore, we have shown that 3bSB can be superior to 3SA. These results indicate that higher-order SB can be useful for solving combinatorial optimization problems expressed naturally by higher-order cost functions. Since we have formulated the higher-order SB by matrix-vector operations, we expect that the higher-order SB may be accelerated by parallel processing with, e.g., GPUs and FPGAs.

We thank Y. Sakai, T. Kashimata, K. Tatsumura, and R. Hidaka for helpful discussions.