Circumventing spin glass traps by microcanonical spontaneous symmetry breaking

We first implement a simple MMC algorithm to explore configurations at the vicinity of an objective energy $E_{\rm o}\!\equiv\!Nu$, whose energy density $u\!\in\!(u_{\rm d},u_{\rm mic})$. An initial configuration $\vec{\sigma}$ with energy $E(\vec{\sigma})\!\leq\!E_{\rm o}$ can be sampled according to the equilibrium Boltzmann distribution (6) at a suitably chosen inverse temperature $\beta$. Starting from this initial configuration, the following single-flip trial is performed at each elementary time interval $1/N$ of the MMC evolution: (1) a vertex $j$ is picked uniformly at random from the $N$ vertices and the flip $\sigma_{j}\!\rightarrow\!-\sigma_{j}$ is proposed; (2) if the new configuration energy does not exceed $E_{\rm o}$, this flip is accepted, otherwise it is rejected and the old spin $\sigma_{j}$ is kept for vertex $j$. One sweep (unit time) of this MMC process corresponds to $N$ consecutive single-flip trials.

This single-flip dynamics obeys detailed balance since the forward and backward transitions between two configurations are equally likely Creutz (1983); Rose and Machta (2019). This microcanonical detailed balance condition ensures that all the visited configurations along the evolution trajectory share the same statistical weight. We sample configurations at fixed time intervals and record their overlap values $m$ and energies $E$. The microcanonical inverse temperature of the system at energy density $u$ is estimated through

where $\bigl{\langle}E(\vec{\sigma})\bigr{\rangle}$ is the averaged energy of the sampled configurations Creutz (1983).

Our MMC dynamics is an unbiased random walk within the microcanonical configuration space. The evolution trajectory is expected to be ergodic since the energy density $u$ is above the dynamical spin glass transition point $u_{\rm d}$. As evolution time goes to infinity every configuration of this space should have equal frequency to be visited, and the system will then certainly be in the MP phase since this phase is exponentially dominant in statistical weight. Empirically, we have observed that this MMC dynamics does indeed achieve MSSB in small-sized systems [Fig. 2 and Fig. 3], while the evolution time needed to witness such a transition increases quickly with system size $N$ [Fig. 3]. It is easier for the MSSB transition to occur in networks of larger degrees $K$. This is consistent with the mean-field theoretical results, which reveal that the entropy barrier and the gap of the order parameter $m$ both decreases with $K$ (Table 1). The hardest RR problem instances are those with degree $K\!=\!4$.

The MSSB transition is an entropy-barrier crossing process [Fig. 3, and the route $A\!\rightarrow\!D\!\rightarrow\!E$ in Fig. 4]. We can regard the MMC dynamics on the coarse-grained level of overlap $m$ as a one-dimensional random walk under a potential field $-s_{u}(m)$, where $s_{u}(m)$ is the entropy density of configurations at fixed energy density $u$ having overlap $m$ with the planted configuration $\vec{\sigma}^{0}$. At each elementary trial, $m$ may change to $m\pm 2/N$ with probability proportional to $(1\mp g_{m})/2$. Here the bias ratio is

with $s_{u}^{\prime}(m)\!\equiv\!{\rm d}s_{u}(m)/{\rm d}m$ being the slope of $s_{u}(m)$ at $m$. Because $s^{\prime}(m)<0$ in $m\!\in\!(0,m^{*})$ with $m^{*}$ being the watershed point of $s_{u}(m)$ [Fig. 3, and point $D$ of Fig. 4], the occurrence of MSSB is an exponentially rare first-passage event characterized by a waiting time of exponential order $O(e^{N\Delta s})$, where $\Delta s\!=\!s_{u}(0)-s_{u}(m^{*})$ is the barrier height Weiss (1981); Khantha and Balakrishnan (1983). Consequently, this simple random walk strategy is practical only for systems of sizes

The entropy barrier $\Delta s$ decreases with energy density $u$ [Fig. 4]. Therefore, we can decrease the entropy barrier $\Delta s$ by fixing $u$ of the MMC dynamics to a lower value. This will help reducing the waiting time of the MSSB transition. On the other hand, at lower energy density values the number of transition paths from point $A$ to point $D$ of Fig. 4 becomes much more rarer and then disappears at $u\!\approx\!u_{\rm d}$. This later kinetic effect will lead to an increase in the waiting time. There should be an optimal energy density $u$, at which the MSSB transition is easiest to observe. We leave this interesting issue to future studies.

As a minor point, we notice that the watershed entropy density $s_{u}(m^{*})$ of Fig. 3 decreases as the energy density $u$ decreases and it will reach zero at certain critical energy density value (denoted as $u_{\rm meb}$) which marks the point of strong microcanonical ergodicity-breaking. When $u\!<\!u_{\rm meb}$ there will be an extended region of the intermediate overlap values $m$ within which there is no any configuration of energy density $u$. In other words, it is impossible to reach from a disordered initial configuration to a highly ordered final configuration through a sequence of small changes of the configuration at fixed energy density $u\!<\!u_{\rm meb}$. We find that $u_{\rm med}$ is below $u_{\rm d}$ for all the degrees $K\!\geq\!3$ (see Figs. 2, 7 and 8), so this strong microcanonical ergodicity-breaking is not really relevant for our MSSB mechanism which occurs at $u\!>\!u_{\rm d}$.

Besides the microcanonical route $A\!\rightarrow\!D\!\rightarrow\!E$ as depicted in Fig. 4, there is also the conventional canonical route ($A\!\rightarrow\!B\!\rightarrow\!C$) to escape the disordered symmetric phase. Will this canonical route be blocked by the low-energy disordered spin glass states? If the planted ground state could be successfully approached following either the microcanonical route or the canonical route, which of these two routes will be more efficient? We explore these two questions in this and the next subsections.

To be fair in comparing the MMC and CMC routes, we design a hybrid-flip Monte Carlo dynamics which is rejection-free both at fixed energy density $u$ (MMC) and at fixed inverse temperature $\beta$ (CMC). Besides the action of flipping a single spin, the hybrid-flip dynamics differs from the single-flip dynamics by allowing the action of flipping a pair of spins in one elementary updating step.

An elementary step of the hybrid-flip MMC dynamics goes as follows: (1) a vertex $j$ is picked uniformly at random from the $N$ vertices and its spin $\sigma_{j}$ is flipped, causing a change $\Delta E$ to the energy; (2) if and only if the energy of the new configuration exceeds the objective value $E_{\rm o}$, an additional flip is performed on the spin $\sigma_{k}$ of a vertex $k$ to bring the energy back to its original value ($k$ is picked uniformly at random from all the vertices capable of causing a change $-\Delta E$ to the energy). This microcanonical hybrid-flip rule obeys microcanonical detailed balance, so all the visited configurations on the MMC evolution trajectory have the same statistical weight.

An elementary step of the hybrid-flip CMC dynamics is quite similar: (1) a vertex $j$ is picked uniformly at random from the $N$ vertices and its spin $\sigma_{j}$ is flipped, causing a change $\Delta E$ to the energy; (2) if $\Delta E$ is positive, then with probability $1\!-\!e^{-\beta\Delta E}$ an additional flip is performed on the spin $\sigma_{k}$ of a vertex $k$ to bring the energy back to its original value ($k$ is picked uniformly at random from all the vertices capable of causing a change $-\Delta E$ to the energy). This canonical hybrid-flip rule obeys canonical detailed balance and the visited configurations on the CMC evolution trajectory are governed by the Boltzmann distribution (6).

We carry out hybrid-flip MMC and hybrid-flip CMC simulations on a single random graph instance of degree $K\!=\!10$ and size $N\!=\!1008$. A total number of $4000$ independent MMC evolution trajectories and $4000$ independent CMC evolution trajectories are sampled with the help of the TRNG pseudo-random number library Bauke and Mertens (2007), all of these trajectories start from the same initial equilibrium configuration which is located in the DS phase (inverse temperature $\beta\!=\!0.4847$ and energy density $u\!=\!-1.50$). We regard the DS phase as successfully escaped when the overlap $m$ with the planted configuration exceeds $0.6$ for the first time, and the evolution is then terminated and the total number $W$ of spin flips along the whole evolution trajectory is returned. We refer to $W$ as the effort to escape DS. Figure 4 shows the cumulative distribution profile $C(W)$ of the effort $W$ over the $4000$ samples, obtained at fixed inverse temperature $\beta$ (CMC) or fixed energy density $u$ (MMC).

We find that the effort (number of spin flips) $W$ fluctuates over a wide range from $10^{6}$ up to $10^{12}$ among repeated runs (all of which start from the same initial configuration). This indicates that restart might be an optimal strategy if the evolution trajectory fails to escape the DS phase within a certain threshold value of $W$. The cumulative distribution function $C(W)$ is concave in shape, both for the MMC dynamics and the CMC dynamics. For example, within an effort $W_{1}\!\approx\!6.071\times 10^{9}$ the success probability of escaping the DS phase is $C(W_{1})\!=\!20\%$ by the CMC dynamics, while this success probability only increases to $C(W_{2})\!=\!35.6\%$ when the effort is doubled to $W_{2}\!=\!1.214\times 10^{10}$ [Fig. 4]. This concave property of $C(W)$ guarantees that restart is an optimal strategy. For example, we could restart the CMC or the MMC evolution from the fixed initial configuration if the effort of the current trajectory exceeds $W_{1}$, and the resulting cumulative distribution of the effort will change to be approximately linear.

Our CMC simulation results clearly demonstrate that the low-energy spin glass states are not encountered by the canonical route $A\!\rightarrow\!B\!\rightarrow\!C$. This may be easy to understand: the canonical route occurs at a much lower inverse temperature than the critical value $\beta_{\rm d}$ where the spin glass dynamical phase transition occurs, and the energy density of the watershed point $B$ separating the DS and the CP phases is considerably higher than the critical value $u_{\rm d}$ of spin glass dynamical phase transition.

Figure 4 demonstrates that, when starting from the same initial configuration, the CMC dynamics is more efficient than the MMC dynamics. For example, the MMC dynamics needs a larger effort $W_{1}\!\approx\!7.093\times 10^{9}$ to achieve a success probability $C(W_{1})\!=\!20\%$. The empirical observation of CMC being more efficient than MMC is consistent with the theoretical predictions on the entropy-barrier height $\Delta s$ and the rescaled free energy barrier height $\beta\Delta f$ [Fig. 4]. We could prove that $\Delta s>\beta\Delta f$ as follows.

The rescaled free energy barrier ($\beta\Delta f$) of the canonical route $A\!\rightarrow\!B\!\rightarrow\!C$ is

Here $\beta_{X}$, $u_{X}$, $f_{X}$, and $s_{X}$ denote the inverse temperature, energy density, free energy density, and entropy density at the phase point $X\in\{A,B,\ldots\}$ of Fig. 4. The entropy barrier ($\Delta s$) of the microcanonical route $A\!\rightarrow\!D\!\rightarrow\!E$) is

Becuase of the relationship (8) and the fact that $\beta_{B}=\beta_{A}$ and $u_{D}=u_{A}$, we know that

In deriving this inequality, we have used the following approximation

assuming that $u_{D}$ is quite close to $u_{B}$.

On the other hand, if we compare the CMC and MMC routes passing through the same intermediate state (say point $B$ of Fig. 4), namely the CMC route $A\!\rightarrow\!B\!\rightarrow\!C$ and the MMC route $A^{\prime}\!\rightarrow B\!\rightarrow\!E^{\prime}$ (points $A^{\prime},B,E^{\prime}$ have the same energy density), we find that the entropy barrier

is lower than the rescaled free energy barrier $\beta_{B}\Delta f$ [Eq. (12)]. To prove this, we notice that

as $\beta_{A^{\prime}}\!>\!\beta_{A}$ and $u_{A}\!>\!u_{A^{\prime}}$. This indicates that for all the transition trajectories passing through the same intermediate state $B$, the MMC trajectories are more efficient than the CMC trajectories. Our simulation results confirm this prediction, see the dotted line and the solid line of Fig. 4.

We do not need to fix the energy or the inverse temperature in the Monte Carlo simulations. Instead the energy or the inverse temperature can slowly change with time. Following some previous numerical efforts Alava et al. (2008); Ma and Zhou (2011); Rose and Machta (2019), we implement SEA (simulated energy annealing) as follows: (1) Start from an initial configuration of high energy value $E$; (2) perform MMC simulation with at objective energy $E_{\rm o}=E$ for a total number $\tau$ of sweeps (each sweep corresponds to $N$ contiguous hybrid-flip updates); (2) decrease the objective energy by a small amount (say, $E_{\rm o}\!-\!2\rightarrow E_{\rm o}$) and repeat step (2).

Some SEA evolution trajectories are shown in Fig. 5 in the $(u,\beta)$ and $(u,m)$ phase planes. As the waiting time $\tau$ at each energy density $u$ increases the MSSB transition occurs at higher energy values $u$. This is consistent with expectations.

Similar results are also observed in the conventional simulated temperature annealing process, where the canonical inverse temperature $\beta$ gradually increases (but does not exceed $\beta_{\rm d}$). The energy density $u$ will change dramatically during this annealing process.

In comparison with the SEA algorithm, the MMC algorithm (with restart) at a suitably chosen fixed energy density $u$ is more efficient.

Since it is exceedingly difficult to escape the DS phase by both the MMC dynamics and the CMC dynamics as the system size $N$ increases, might it be possible to retrieve the planted ground state without actually leaving the DS phase? Here we report a preliminary exploration of this interesting question.

Given a $p$-spin problem instance, it is straightforward to sample a large set of independent equilibrium configurations of the DS phase at a given energy density $u\!\in\!(u_{\rm d},u_{\rm mic})$ by the hybrid-flip MMC algorithm. The overlaps $m$ of these configurations are of order $1/\sqrt{N}$, so a more convenient quantity is the rescaled overlap $\lambda\equiv m\sqrt{N}$, namely

An empirical probability distribution $P(\lambda)$ could be constructed as

where $\mathcal{N}$ denotes the total number of samples, $\vec{\sigma}^{(\ell)}=(\sigma_{1}^{(\ell)},\ldots,\sigma_{N}^{(\ell)})$ is the $\ell$-th sample, and $\delta(x)\!=\!1$ for $x\!=\!0$ and $\delta(x)\!=\!0$ for $x\!\neq\!0$.

For the planted $3$-spin systems of degree $K\!\geq\!4$ at zero noise ($\varepsilon\!=\!0$), we find that the overlap distribution $P(\lambda)$ is not symmetric in $\lambda$, see Fig. 6 for some representative examples. When the magnitude $|\lambda|$ is small the rescaled overlap $\lambda$ is more likely to be negative than to be positive, while it is more likely to be positive than to be negative as $|\lambda|$ increases beyond certain moderate value [Fig. 6]. Overall, the rescaled overlap $\lambda$ is more likely to be negative than to be positive. Our extensive empirical results suggest that the probabilities of $\lambda$ being non-positive ($P_{\leq 0}$) and non-negative ($P_{\geq 0}$) scale with $N$ as

where the asymmetry parameter $\gamma$ is estimated to be $\gamma\approx 1.0$ through the bootstrap method Efron (1979). (Very surprisingly, our empirical results on the first moment $\langle\lambda\rangle$ suggest that it is identical to zero.)

We find that the empirical $P(\lambda)$ distribution could be precisely fitted by

with $U(\lambda)$ being an odd function of the following form

Both the linear-order coefficient $c_{1}$ and the third-order coefficient $c_{3}$ are positive [Fig. 6]. According to the pioneering theoretical work in the 1990s Van den Broeck and Reimann (1996); Reimann and Van den Broeck (1996), the anti-symmetric “potential energy” $U(\lambda)$ may make it possible to infer the planted ground state $\vec{\sigma}^{0}$ by unsupervised learning.

We can formulate an unsupervised learning task for the planted $p$-spin model as follows: (1) Assume all the observed configurations $\vec{\sigma}^{(\ell)}$ are generated by the following generative model

with $N$ Ising parameters $\{\sigma_{i}^{0}\}$. (2) Determine the optimal assignment for this set of model parameters so that these $\mathcal{N}$ observed configurations will be most likely to be observed Engel and Van den Broeck (2001).

We expect the planted ground state to be the true optimal solution of this unsupervised learning problem. But it is also likely that this problem has many sub-optimal solutions, some of which being the spin glass minimal-energy configurations of Eq. (1). Whether the optimal solution could be efficiently achieved for this machine-learning problem is an open issue and deserves thorough investigations in the future.

Other statistical inference strategies, such as perceptron learning and principal components analysis may also help for the planted $p$-spin model Engel and Van den Broeck (2001); Kabashima and Uda (2004); Braunstein and Zecchina (2006); Huang and Toyoizumi (2015); Zhou (2019b); Hu et al. (2019). For example, the asymmetric property (20) of the sampled configurations indicates the existence of an optimal hyperplane separating the $\mathcal{N}$ configuration samples in the most asymmetric manner.