Combinatorial optimization using dynamical phase transitions in driven-dissipative systems

The dynamics of driven-dissipative systems is shown to be well-fitted for achieving efficient combinatorial optimization. The proposed method can be applied to solve any combinatorial optimization problem that is equivalent to minimizing an Ising Hamiltonian. Moreover, the dynamics considered can be implemented using various physical systems as it is based on generic dynamics—the normal form of the supercritical pitchfork bifurcation. The computational principle of the proposed method relies on an hybrid analog-digital representation of the binary Ising spins by considering the gradient descent of a Lyapunov function that is the sum of an analog Ising Hamiltonian and archetypal single or double-well potentials. By gradually changing the shape of the latter potentials from a single to double well shape, it can be shown that the first nonzero steady states to become stable are associated with global minima of the Ising Hamiltonian, under the approximation that all analog spins have the same amplitude. In the more general case, the heterogeneity in amplitude between analog spins induces the stabilization of local minima, which reduces the quality of solutions to combinatorial optimization problems. However, we show that the heterogeneity in amplitude can be reduced by setting the parameters of the driving signal near a regime, called the dynamic phase transition, where the analog spins' DC components map more accurately the global minima of the Ising Hamiltonian which, in turn, increases the quality of solutions found. Last, we discuss the possibility of a physical implementation of the proposed method using networks of degenerate optical parametric oscillators.

Figure 1

(a) Landscape of a configuration-dependent function $\frac{2\epsilon }{N}\mathcal{H}$ and the bifurcation parameter $\alpha$ in the case of homogeneous amplitudes, i.e., $x_j^2=x^2, \forall j$. (b) The same as (a) (see dotted line), but the case of heterogeneous amplitudes (see full line) is superimposed. Note that the first nonzero stable steady-state configuration does not always correspond to the ground-state configuration of the Ising Hamiltonian in this case.

Figure 2

(a) Driven-dissipative system in the symmetry-restoring oscillation (SRO) and symmetry-breaking oscillation (SBO) phases shown by the solid and dashed lines, respectively. (b) Double-well potential $V_b$. Upper and lower arrows denote a schematic representation of SRO and SBO, respectively.

Figure 3

(a) Average absolute amplitude $\langle \left|\overline{x}\right|\rangle$ of the signals ${\left\{x_i\right\}}_i$ shown by the color gradient with $\langle |\overline{x}|\rangle =\frac{1}{N}{\sum }_j\left|\frac{1}{T}{\int }_0^T{x}_j\left(t\right)dt\right|$ and $T$ the duration of the simulation shown in the parameter space $(h,f)$. The system is driven by a random phase shift key modulation [see Eq. (16)]. The solid and dotted lines show the set of parameters at the dynamic phase transition, i.e., for which ${\tilde{\lambda }}_1=0$, in the coupled and uncoupled case, respectively. (b) The same as (a), but for the minimal value of the Ising Hamiltonian $\mathcal{H}$ calculated using the instantaneous values of the spin ${\sigma }_j=\frac{x_j}{|x_j|}$. $N=64, \epsilon =0.5$, and $p=8$.

Figure 4

Average trajectory of the Ising Hamiltonian $\mathcal{H}$ shown by solid and dashed curves in the case of a random phase shift keying modulation case and a nondriven case, respectively. The minimum value of $\frac{\mathcal{H}}{N}$ is calculated over each period of the driving signal in the three-dimensional Edwards-Anderson model with $N=216$ and normally distributed interactions. $\epsilon =0.5, h=100$, and $f\approx 21.7$. $\alpha =1$ and $\alpha =-0.3$ in the driven and nondriven case, respectively.

Figure 5

Reduction of the amplitude heterogeneity due to the driving signal. For the combinatorial optimization problem considered in this figure, the ground states and first excited states have Ising Hamiltonian equal to $\mathcal{H}=-8$ and $\mathcal{H}=-6$, respectively. (a1, a2) Network structure $w_{ij}$, spin directions, and analog spins' squared amplitudes are shown by lines, arrows, and circles, respectively. The radii of circles are proportional to the squared amplitudes $x_j^2$ calculated using a single run. The combinatorial optimization problem encoded by the network structure shown in (a1) and (a2) is such that the connections between spins are $\epsilon w_{ij}$ with $w_{ij}\in \{-1,0\}$ and $\epsilon =0.1$. In (a1) and (a2) are shown the nondriven and driven cases, respectively, for the bifurcation parameter $\alpha$ such that $\alpha -{\alpha }_c=0.1$, where ${\alpha }_c$ is the critical value at the corresponding phase transitions. (b1) Probability $P_0$ to find the ground-state configurations vs. $\alpha -{\alpha }_c$. (b2) Averaged amplitude heterogeneity ${\delta }_{\overline{x}}$ vs. $\alpha -{\alpha }_c$. In (b1, b2), the nondriven and driven cases are shown by full and dashed curves, respectively, and 200 runs are used for the calculation of $P_0$ and ${\delta }_{\overline{x}}$.

Figure 6

(a) Average absolute amplitude $\langle \left|\overline{x}\right|\rangle$ (see definition in caption of Fig. 3) and eigenvalue ${\tilde{\lambda }}_1$ in the case of a random phase shift keying signal shown by the thick and thin curves, respectively. The solid, dotted, and dashed curves show the cases $c_{ij}=0, c_{ij}\propto {\omega }_{ij}$, and $c_{ij}\propto -{\omega }_{ij}$, respectively. (b) Amplitude heterogeneity ${\delta }_x$ vs. modulation frequency $f-f_c$ with $f_c$ the frequency at the dynamic phase transition [see Eq. (38)]. (c) The minimal value of the Ising Hamiltonian $\text{min}\left(\frac{\mathcal{H}}{N}\right)$ vs. modulation frequency $f-f_c$ in the three-dimensional Edwards-Anderson model with $N=216$ and normally distributed interactions. $\epsilon =0.5, \alpha =7$, and $h=100$.

Figure 7

(a) The minimal Ising Hamiltonian $\text{min}\left(\frac{\mathcal{H}}{N}\right)$ vs. the bifurcation parameter $\alpha$ at the optimal driving frequency $f^{*}$ in the case of uncorrelated and correlated [see Eq. (41)] driving signals shown by the solid and dotted curves, respectively. Bars show the $95\%$ confidence intervals. $f^{*}\approx 7.7, f^{*}\approx 10.2, f^{*}\approx 12.4$, and $f^{*}\approx 22$ for $\alpha =17, \alpha =11, \alpha =7$, and $\alpha =1$, respectively. $\epsilon =0.5, h=100$, and $T=6\times {10}^2$. (b) The minimal Ising Hamiltonian $min\left(\frac{\mathcal{H}}{N}\right)$ vs. the amplitude $h$ at the optimal driving frequency. $f^{*}\approx 4.55, f^{*}\approx 10.2, f^{*}\approx 15.6$, and $f^{*}\approx 20.8$ for $h=50, h=100, h=150$, and $h=200$, respectively. $\epsilon =0.5, \alpha =11$, and $T=3\times {10}^2$.

Figure 8

(a, b) Minimal values of the Ising Hamiltonian $\mathcal{H}$ as a function of the problem size. In (a) and (b) are shown the cases of normally distributed and $\pm 1$ interactions, respectively. $\epsilon =0.5, h=100$, and $f\approx 12$. In the driven and nondriven cases, the pump rate is $\alpha =7$ and $\alpha \approx -0.8$, respectively. (a) $T=3\times {10}^2$. (b) $T=8\times {10}^2$.

Figure 9

The minimum value of the Ising Hamiltonian $\frac{\mathcal{H}}{N}$ vs. the driving frequency $f^{*}$ in the driven cases $c_{ij}=0, c_{ij}\propto {\omega }_{ij}$, and $c_{ij}=\pm 1$ shown in (a), (b), and (d), respectively. $h=100, \alpha =7$, and $\epsilon =0.5$. (c) The minimum value of the Ising Hamiltonian $\frac{\mathcal{H}}{N}$ vs. the bifurcation parameter $\alpha$ in the nondriven case $d_i=0$. $\epsilon =0.5$.

Figure 10

Performance of the proposed method in finding the ground state of the Edwards-Anderson model for each problem size $N$ shown in the uncorrelated, correlated, and nondriven cases by dots, diamonds, and squares, respectively. (a1, a2) Probability $P_0$ to find the ground state for a duration of $T$. (b1, b2) Averaged computational time $t_1$ to the ground state for successful runs. (c1, c2) Averaged computational time $t_0$ required to find the ground state with 99% success probability. $\epsilon =0.5, \alpha =17, h=100, f\approx 8$. In (a1, b1, c1) and (a2, b2, c2) are shown the cases of normally distributed and $\pm 1$ interactions, respectively. The lines superimposed to the data are obtained using linear regression in the linear-logarithmic scale. (a1, b1, c1) $T=8\times {10}^2$ for $N=64, N=125$, and $N=216$. $T=16\times {10}^2$ for $N=343$. (a2, b2, c2) $T=8\times {10}^2$.

Figure 11

(a, b) Probability $P_0$ of finding the ground state vs. the system size $N$ in the driven and nondriven cases, respectively, for several parameters $\epsilon$. $T={10}^3, \alpha =0.5, f^{*}=0.1$, and weights given by the three-dimensional EA model with ${\omega }_{jl}=\pm 1$.