Efficiency of quantum versus classical annealing in non-convex learning problems

TL;DR

This paper compares quantum and classical annealing methods for non-convex learning problems, showing quantum annealing's efficiency in escaping local minima where classical methods struggle, especially in machine learning contexts.

Contribution

It identifies classes of non-convex problems where quantum annealing outperforms classical thermal annealing, highlighting its potential in machine learning optimization tasks.

Findings

Quantum annealing efficiently finds solutions in complex landscapes.

Classical thermal annealing experiences exponential slowdowns.

Quantum tunneling enables escape from local minima.

Abstract

Quantum annealers aim at solving non-convex optimization problems by exploiting cooperative tunneling effects to escape local minima. The underlying idea consists in designing a classical energy function whose ground states are the sought optimal solutions of the original optimization problem and add a controllable quantum transverse field to generate tunneling processes. A key challenge is to identify classes of non-convex optimization problems for which quantum annealing remains efficient while thermal annealing fails. We show that this happens for a wide class of problems which are central to machine learning. Their energy landscapes is dominated by local minima that cause exponential slow down of classical thermal annealers while simulated quantum annealing converges efficiently to rare dense regions of optimal solutions.

Tables

Table 1. Table 1: Results for the 20 20 20 small samples (original an randomized) at the end of the QA process. The samples are the same as for Figg.  8 and  9 , where they are arranged in row-major order. The second column shows the number of solutions; the other columns are described in the text.

sample	#sol	Original samples
		$P_{\mathrm{SOL}}$	$⟨E⟩/N$	$E\left({\sigma }^{\star }\right)$	$p\left({\sigma }^{\star }\right)$	$\mathrm{IPR}$	$\overline{d}$	$r_{\mathrm{SQA}}$
1	80	0.997873	0.007283	0	0.091049	0.052860	0.161271	2
2	152	0.999345	0.005566	0	0.077510	0.035271	0.171339	30
3	25	0.998965	0.009524	0	0.260645	0.132882	0.116819	5
4	129	0.999449	0.005236	0	0.163546	0.067067	0.121317	18
5	41	0.989324	0.022258	0	0.295728	0.149694	0.109064	1
6	24	0.964274	0.041914	0	0.287578	0.169719	0.120534	1
7	24	0.998103	0.010600	0	0.304858	0.150627	0.150411	3
8	40	0.999237	0.006526	0	0.165771	0.078460	0.168933	3
9	48	0.999253	0.005652	0	0.146795	0.070516	0.122878	2
10	149	0.999120	0.009842	0	0.166722	0.084118	0.129071	16
11	27	0.999517	0.003968	0	0.292468	0.161192	0.101294	1
12	69	0.999140	0.009151	0	0.325233	0.151204	0.081271	12
13	47	0.999550	0.004662	0	0.424602	0.230421	0.062514	3
14	56	0.999123	0.006101	0	0.316729	0.151897	0.136677	8
15	54	0.999430	0.005455	0	0.192050	0.100316	0.132264	1
16	28	0.994853	0.010547	0	0.175557	0.115669	0.184412	3
17	49	0.999546	0.004504	0	0.344131	0.169363	0.122593	5
18	28	0.999361	0.006120	0	0.359053	0.187311	0.103266	6
19	41	0.998693	0.004844	0	0.244020	0.111271	0.136971	1
20	22	0.997396	0.007165	0	0.161562	0.108546	0.098991	4
mean	56.65	0.996578	0.009346	0.0	0.239781	0.123921	0.126595	-

Equations

$$E\left(\left\{{\sigma }_j\right\}\right)=\sum _{\mu =1}^{\alpha N}{\Delta }_{\mu }^n\Theta \left(-{\Delta }_{\mu }\right),{\Delta }_{\mu }\doteq \frac{{\tau }^{\mu }}{\sqrt{N}}\sum _{j=1}^N{\xi }_j^{\mu }{\sigma }_j$$

$$\hat{H}=E\left(\left\{{\hat{\sigma }}_j^z\right\}\right)-\Gamma \sum _{j=1}^N{\hat{\sigma }}_j^x$$

$$H_{\mathrm{eff}}\left({\left\{{\sigma }_j^a\right\}}_{j,a}\right)=\frac{1}{y}\sum _{a=1}^{y}E\left({\left\{{\sigma }_j^a\right\}}_j\right)-\frac{\gamma }{\beta }\sum _{a=1}^y\sum _{j=1}^N{\sigma }_j^a{\sigma }_j^{a+1}-\frac{NK}{\beta }$$

$$H_{\mathrm{eff}}^{\mathrm{RE}}\left({\sigma }^{\star },{\left\{{\sigma }_j^a\right\}}_{j,a}\right)=\sum _{a=1}^{y}E\left({\left\{{\sigma }_j^a\right\}}_j\right)-\frac{\lambda }{\beta }\sum _{a=1}^y\sum _{j=1}^N{\sigma }_j^a{\sigma }_j^{\star }$$

$$\hat{H}=E\left({\left\{{\hat{\sigma }}_j^z\right\}}_j\right)-\Gamma \sum _{j=1}^N{\hat{\sigma }}_j^x$$

$$Z=\mathrm{Tr}\left(e^{-\beta \hat{H}}\right).$$

$$H_{\mathrm{eff}}\left({\left\{{\sigma }_j^a\right\}}_{j,a}\right)=\frac{1}{y}\sum _{a}E\left({\left\{{\sigma }_j^a\right\}}_j\right)-\frac{\gamma }{\beta }\sum _{aj}{\sigma }_j^a{\sigma }_j^{a+1}-\frac{NK}{\beta }$$

$$\gamma$$

$$K$$

$$Z_{\mathrm{eff}}=\sum _{\left\{{\sigma }^a\right\}}{e}^{-\frac{\beta }{y}{\sum }_{a}E\left({\sigma }^a\right)+\gamma {\sum }_{aj}{\sigma }_j^a{\sigma }_j^{a+1}+NK}.$$

$$E\left(\sigma \right)=\sum _{\mu =1}^{\alpha N}\Theta \left(-\frac{1}{\sqrt{N}}\sum _j{\xi }_j^{\mu }{\sigma }_j\right)$$

$${⟨\log Z⟩}_{\xi }$$

$${⟨Z_{\mathrm{eff}}^n⟩}_{\xi }=e^{nNK}{⟨\int \prod _{caj}d\mu \left({\sigma }_j^{ca}\right)\prod _{caj}{e}^{\gamma {\sigma }_j^c{\sigma }_j^{c\left(a+1\right)}}\prod _{\mu ca}\left(\Theta \left[\frac{1}{\sqrt{N}}\sum _j{\xi }_j^{\mu }{\sigma }_j^{ca}\right]\left(1-e^{-\frac{\beta }{y}}\right)+e^{-\frac{\beta }{y}}\right)⟩}_{\xi }$$

$${⟨Z_{\mathrm{eff}}^n⟩}_{\xi }$$

$${⟨\prod _{\mu ca}\delta \left({\lambda }_{\mu }^{ca}-\frac{1}{\sqrt{N}}\sum _j{\xi }_j^{\mu }{\sigma }_j^{ca}\right)⟩}_{\xi }=$$

$$\int \prod _{\mu ca}\frac{d{\hat{\lambda }}_{\mu }^{ca}}{2\pi }\prod _{\mu ca}{e}^{i{\hat{\lambda }}_{\mu }^{ca}{\lambda }_{\mu }^{ca}}\prod _{\mu }\exp \left(-\frac{1}{2}\sum _{cdab}{\hat{\lambda }}_{\mu }^{ca}{\hat{\lambda }}_{\mu }^{db}\left(\frac{1}{N}\sum _j{\sigma }_j^{ca}{\sigma }_j^{db}\right)\right)$$

$${⟨Z_{\mathrm{eff}}^n⟩}_{\xi }$$

$$G_S$$

$$G_E$$

$$q^{ca,db}=\{\begin{array}{ll}{q}_1& \mathrm{if}c=d\\ q_0& \mathrm{if}c\ne d\end{array}$$

$${⟨Z_{\mathrm{eff}}^n⟩}_{\xi }$$

$$G_S$$

$$G_E$$

$$G_S$$

$$g\left(z_{0,}{z}_1,w\right)$$

$$h\left(z_0,z_1\right)$$

$${\mathcal{G}}_S$$

$$G_E$$

$${\mathcal{G}}_E\doteq \frac{1}{n}\log G_E=\int D{z}_0\log \int D{z}_1\exp \left(-\beta H\left(\frac{z_1\sqrt{q_1-q_0}+z_0\sqrt{q_0}}{\sqrt{1-q_1}}\right)\right)$$

$$\varphi \doteq \frac{1}{N}⟨\log Z_{\mathrm{eff}}⟩={\mathrm{extr}}_{q_0,q_1,{\hat{q}}_0,{\hat{q}}_1}\left\{\frac{1}{2}{y}^2{q}_0{\hat{q}}_0-\frac{1}{2}y\left(y-1\right)q_1{\hat{q}}_1-\frac{1}{2}{\hat{q}}_{1}y+{\mathcal{G}}_S+\alpha {\mathcal{G}}_E\right\}$$

$${\hat{q}}_0$$

$${\hat{q}}_1$$

$$\varphi$$

$${\mathcal{G}}_S$$

Full text

1. I Introduction

2. II Energy functions

3. III Connection with the local entropy measure

4. IV Phase diagram: analytical and numerical results

5. V Conclusions

6. A Theoretical analysis by the replica method

7. A.0.1 Small $\Gamma$ limit

8. A.1 Energy function with stability

9. A.1.1 Small $\Gamma$ limit

10. B Estimation of the local energy and entropy landscapes with the cavity method

11. C Numerical simulations details of the annealing protocols

12. C.1 Quantum annealing protocol

13. C.2 Classical simulated annealing protocol

14. D Additional numerical results on the annealing processes

15. D.1 Additional comparisons between theory and simulations

16. D.2 Experiments with two-layer networks

17. E Real-time Quantum Annealing on small samples

18. E.1 Numerical methods

19. E.2 Sample selection

20. E.2.1 Randomized samples

21. E.3 Analysis

22. E.3.1 QA vs SQA

23. E.3.2 Other measurements

24. E.3.3 Local entropies

25. E.3.4 Energy gaps

I Introduction

Quantum tunneling and quantum correlations govern the behavior of very complex collective phenomena in quantum physics at low temperature. Since the discovery of the factoring quantum algorithms in the 90s Shor (1994), a lot of efforts have been devoted to the understanding of how quantum fluctuations could be exploited to find low-energy configurations of energy functions which encode the solutions of non-convex optimization problems in their ground states. This has led to the notion of controlled quantum adiabatic evolution, where a time dependent many-body quantum system is evolved towards its ground states so as to escape local minima through multiple tunneling events Ray et al. (1989); Finnila et al. (1994); Kadowaki and Nishimori (1998); Farhi et al. (2001); Das and Chakrabarti (2008). When finite temperature effects have to be taken into account, the computational process is called Quantum Annealing (QA). Classical Simulated Annealing (SA) uses thermal fluctuations for the same computational purpose, and Markov Chains based on this principle are among the most widespread optimization techniques across science Moore and Mertens (2011). Quantum fluctuations are qualitatively different from thermal fluctuations and in principle quantum annealing algorithms could lead to extremely powerful alternative computational devices.

In the quantum annealing approach, a time dependent quantum transverse field is added to the classical energy function leading to an interpolating Hamiltonian that may take advantage of correlated fluctuations mediated by tunneling. Starting with a high transverse field, the quantum model system can be initialized in its ground state, e.g. all spins aligned in the direction of the field. The adiabatic theorem then ensures that by slowly reducing the transverse field the system remains in the ground state of the interpolating Hamiltonian. At the end of the process the transverse field vanishes and the systems ends up in the sought ground state of the classical energy function. The original optimization problem would then be solved if the overall process could take place in a time bounded by some low degree polynomial in the size of the problem. Unfortunately, the adiabatic process can become extremely slow. The adiabatic theorem requires the rate of change of the Hamiltonian to be smaller than the square of the gap between the ground state and the first excited state Born and Fock (1928); Landau (1932); Zener (1932). For small gaps the process can thus become inefficient. Exponentially small gaps are not only possible in worst case scenarios but have also been found to exist in typical random systems where comparative studies between quantum and classical annealing have so far failed in displaying quantum exponential speed up, e.g. at first order phase transition in quantum spin glasses Altshuler et al. (2010); Bapst et al. (2013) or 2D spin glass systems Santoro et al. (2002); Martoňák et al. (2002); Heim et al. (2015). More positive results have been found for ad hoc energy functions in which global minima are planted in such a way that tunneling cascades can become more efficient than thermal fluctuations Rønnow et al. (2014); Farhi et al. (2001). As far as the physical implementations of quantum annealers is concerned, studies have been focused on discriminating the presence of quantum effects rather than on their computational effectiveness Johnson et al. (2011); Boixo et al. (2014); Langbein et al. (2004).

Consequently, a key open question is to identify classes of relevant optimization problems for which quantum annealing can be shown to be exponentially faster than its classical thermal counterpart.

Here we give an answer to this question by providing analytic and simulation evidence of exponential speed up of quantum versus classical simulated annealing for a representative class of random non-convex optimization problems of basic interest in machine learning. The simplest example of this class is the problem of training binary neural networks (described in detail below): very schematically, the variables of the problem are the (binary) connection weights, while the energy measures the training error over a given dataset.

These problems have been very recently found to possess a rather distinctive geometrical structure of ground states Baldassi et al. (2015, 2016a, 2016b, 2016c): the free energy landscape has been shown to be characterized by the existence of an exponentially large number metastable states and isolated ground states, and a few regions where the ground states are dense. These dense regions, which had previously escaped the equilibrium statistical physics analysis Krauth and Mézard (1989); Sompolinsky et al. (1990), are exponentially rare, but still possess a very high local internal entropy: they are composed of ground states that are surrounded, at extensive but relatively small distances, by exponentially many other ground states. Under these circumstances, classical SA (as any Markov Chain satisfying detailed balance) gets trapped in the metastable states, suffering ergodicity breaking and exponential slowing down toward the low energy configurations. These problems have been considered to be intractable for decades and display deep similarities with disordered spin glass models which are known to never reach equilibrium.

The large deviation analysis that has unveiled the existence of the rare dense regions has led to several novel algorithms, including a Monte Carlo scheme defined over an appropriate objective function Baldassi et al. (2016a) that bears close similarities with a Quantum Monte Carlo (QMC) technique based on the Suzuki-Trotter transformation Das and Chakrabarti (2008). Motivated by this analytical mapping and by the geometrical structure of the dense and degenerate ground states which is expected to favor zero temperature kinetic processes Foini et al. (2010); Biroli and Zamponi (2012), we have conducted a full analytical and numerical statistical physics study of the quantum annealing problem, reaching the conclusion that in the quantum limit the QMC process, i.e. Simulated Quantum Annealing (SQA), can equilibrate efficiently while the classical SA gets stuck in high energy metastable states. These results generalize to multi layered networks.

While it is known that other quasi-optimal classical algorithms for the same problems exist Baldassi et al. (2016a); Hubara et al. (2016); Courbariaux et al. (2015), here we focus on the physical speed up that a quantum annealing approach could provide in finding rare regions of ground states. We provide physical arguments and numerical results supporting the conjecture that the real time quantum annealing dynamics behaves similarly to SQA.

As far as machine learning is concerned, dense regions of low energy configurations (i.e. quasi-flat minima over macroscopic length scales) are of fundamental interest, as they are particularly well-suited for making predictions given the learned data: on the one hand, these regions are by definition robust with respect to fluctuations in a sizable fraction of the weight configurations and as such are less prone to fit the noise. On the other hand, an optimal Bayesian estimate, resulting from a weighted consensus vote on all configurations, would receive a major contribution from one of such regions, compared to a narrow minimum; the centroid of the region (computed according to any reasonable metric which correlates the distance between configurations with the network outcomes) would act as a representative of the region as a whole MacKay (2003). In this respect, it is worth mentioning that in deep learning LeCun et al. (2015) all the learning algorithms which lead to good prediction performance always include effects of a systematically injected noise in the learning phase, a fact that makes the equilibrium Gibbs measure not the stationary measure of the learning protocols and drive the systems towards wide minima. We expect that these results can be generalized to many other classes of non convex optimization problems where local entropy plays a role, ranging from robust optimization to physical disordered systems.

Quantum gate based algorithms for machine learning exist, however the possibility of a physical implementation remains a critical issue Aaronson (2015).

II Energy functions

As a working example, we first consider the problem of learning random patterns in single layer neural network with binary weights, the so called binary perceptron problem Krauth and Mézard (1989). This network maps vectors of $N$ inputs $\xi\in\left\{-1,+1\right\}^{N}$ to binary outputs $\tau=\pm 1$ through the non linear function $\tau=\mathrm{sgn}\left(\sigma\cdot\xi\right)$, where $\sigma\in\left\{-1,+1\right\}^{N}$ is the vector of synaptic weights. Given $\alpha N$ input patterns $\left\{\xi^{\mu}\right\}_{\mu=1}^{\alpha N}$ with $\mu=1,...,\alpha N$ and their corresponding desired outputs $\left\{\tau^{\mu}\right\}_{\mu=1}^{\alpha N}$, the learning problem consists in finding $\sigma$ such that all input patterns are simultaneously classified correctly, i.e. $\mathrm{sgn}\left(\sigma\cdot\xi^{\mu}\right)=\tau^{\mu}$ for all $\mu$. Both the components of the input vectors $\xi_{i}^{\mu}$ and the outputs $\tau^{\mu}$ are independent identically distributed unbiased random variables ($P\left(x\right)=\frac{1}{2}\delta\left(x-1\right)+\frac{1}{2}\delta\left(x+1\right)$). In the binary framework, the procedure for writing a spin Hamiltonian whose ground states are the sought optimal solutions of the original optimization problem is well known Barahona (1982). The energy $E$ of the binary perceptron is proportional to the number of classification errors and can be written as

[TABLE]

where $\Theta\left(x\right)$ is the Heaviside step function: $\Theta\left(x\right)=1$ if $x>0$, $\Theta\left(x\right)=0$ otherwise. When the argument of the $\Theta$ function is positive, the perceptron is implementing the wrong input-output mapping. The exponent $n\in\left\{0,1\right\}$ defines two different forms of the energy functions which have the same zero energy ground states and different structures of local minima. The equilibrium analysis of the binary perceptron problem shows that in the large size limit and for $\alpha<\alpha_{c}\simeq 0.83$ Krauth and Mézard (1989), the energy landscape is dominated by an exponential number of local minima and of zero energy ground states that are typically geometrically isolated Huang and Kabashima (2014), i.e. they have extensive mutual Hamming distances. For both choices of $n$ the problem is computationally hard for SA processes Horner (1992): in the large $N$ limit, a detailed balanced stochastic search process gets stuck in metastable states at energy levels of order $O(N)$ above the ground states.

Following the standard SQA approach, we identify the binary variables $\sigma$ with one of the components of physical quantum spins, say $\sigma^{z}$, and we introduce the Hamiltonian operator of a model of $N$ quantum spins with the perceptron term of Eq. (1) acting in the longitudinal direction $z$ and a magnetic field $\Gamma$ acting in the transverse direction $x$. The interpolating Hamiltonian reads:

[TABLE]

where $\hat{\sigma}_{j}^{z}$ and $\hat{\sigma_{j}}^{x}$ are the spin operators (Pauli matrices) in the $z$ and $x$ directions. For $\Gamma=0$ one recovers the classical optimization problem. The QA procedure consists in initializing the system at large $\beta$ and $\Gamma$, and slowly decreasing $\Gamma$ to [math]. To analyze the low temperature phase diagram of the model we need to study the average of the logarithm of the partition function $Z=\mathrm{Tr}\,\left(e^{-\beta\hat{H}}\right).$ This can be done using the Suzuki-Trotter transformation which leads to the study of a classical effective Hamiltonian acting on a system of $y$ interacting Trotter replicas of the original classical system coupled in an extra dimension:

[TABLE]

where the $\sigma_{j}^{a}=\pm 1$ are Ising spins, $a\in\left\{1,\dots,y\right\}$ is a replica index with periodic boundary conditions $\sigma_{j}^{y+1}\equiv\sigma_{j}^{1}$, $\gamma=\frac{1}{2}\log\coth\left(\frac{\beta\Gamma}{y}\right)$ and $K=\frac{1}{2}y\log\left(\frac{1}{2}\sinh\left(2\frac{\beta\Gamma}{y}\right)\right).$

The replicated system needs to be studied in the limit $y\to\infty$ to recover the so called path integral continuous quantum limit and to make the connection with the behavior of quantum devices Heim et al. (2015). The SQA dynamical process samples configurations from an equilibrium distribution and it is not necessarily equivalent to the real time Schrödinger equation evolution of the system. A particularly dangerous situation occurs if the ground states of the system encounter first order phase transitions which are associated to exponentially small gaps Altshuler et al. (2010); Bapst and Semerjian (2012, 2013) at finite N. As discussed below, this appears not to be the case for the class of models we are considering.

III Connection with the local entropy measure

The effective Hamiltonian Eq. (3) can be interpreted as many replicas of the original systems coupled through one dimensional periodic chains, one for each original spin, see Fig. 1b. Note that the interaction term $\gamma$ diverges as the transverse field $\Gamma$ goes to [math]. This geometrical structure is very similar to that of the Robust Ensemble (RE) formalism Baldassi et al. (2016a), where a probability measure that gives higher weight to rare dense regions of low energy states is introduced. There, the main idea is to maximize $\Phi\left(\sigma^{\star}\right)=\log\sum_{\left\{\sigma\right\}}e^{-\beta E\left(\sigma\right)-\lambda\sum_{j=1}^{N}\sigma_{j}\sigma_{j}^{\star}}$, i.e. a “local free entropy” where $\lambda$ is a Lagrange parameter that controls the extensive size of the region around a reference configuration $\sigma^{\star}$. One can then build a new Gibbs distribution $P\left(\sigma^{\star}\right)\propto e^{y\Phi\left(\sigma^{\star}\right)}$, where $-\Phi$ has the role of an energy and $y$ of an inverse temperature: in the limit of large $y$, this distribution concentrates on the maxima of $\Phi$. Upon restricting the values of $y$ to be integer (and large), $P\left(\sigma^{\star}\right)$ takes a factorized form yielding a replicated probability measure $P_{\mathrm{RE}}\left(\sigma^{\star},\sigma^{1},\dots,\sigma^{y}\right)\propto e^{-\beta H_{\mathrm{eff}}^{\mathrm{RE}}\left(\sigma^{\star},\left\{\sigma_{j}^{a}\right\}\right)}$ where the effective energy is given by

[TABLE]

As in the Suzuki-Trotter formalism, $H_{\mathrm{eff}}^{\mathrm{RE}}\left(\sigma^{\star},\left\{\sigma_{j}^{a}\right\}_{j,a}\right)$ corresponds to a system with an overall energy given by the sum of $y$ individual “real replica energies” plus a geometric coupling term; in this case however the replicas interact with the “reference” configurations $\sigma^{\star}$ rather than among themselves, see Fig. 1c.

The Suzuki-Trotter representation and the RE formalism differ in the topology of the interactions between replicas and in the scaling of the interactions, but for both cases there is a classical limit, $\Gamma\to 0$ and $\lambda\to\infty$ respectively, in which the replicated systems are forced to correlate and eventually coalesce in identical configurations. For non convex problems, these will not in general correspond to configuration dominating the original classical Gibbs measure.

For the sake of clarity we should remind that in the classical limit and for $\alpha<\alpha_{c}$, our model presents an exponential number of far apart isolated ground states which dominate the Gibbs measure. At the same time, there exist rare clusters of ground states with a density close to its maximum possible value (high local entropy) for small but still macroscopic cluster sizes Baldassi et al. (2015). This fact has several consequences: no further subdivision of the clusters into states is possible, the ground states are typically $O(1)$ spin flip connected Baldassi et al. (2015) and a tradeoff between tunneling events and exponential number of destination states within the cluster is possible.

IV Phase diagram: analytical and numerical results

Thanks to the mean field nature of the energetic part of the system, Eq. (3), we can resort to the replica method for calculating analytically the phase diagram. As discussed in the Appendix Sec. A, this can be done under the so called static approximation, which consists in using a single parameter $q_{1}$ to represent the overlaps along the Trotter dimension, $q_{1}^{ab}=\left\langle\frac{1}{N}\sum_{j=1}^{N}\sigma_{j}^{a}\sigma_{j}^{b}\right\rangle\approx q_{1}$. Although this approximation crudely neglects the dependency of $q_{1}^{ab}$ from $\left|a-b\right|$, the resulting predictions show a remarkable agreement with numerical simulations.

In the main panel of Fig. 2, we report the analytical predictions for the average classical component of the energy of the quantum model as a function of the transverse field $\Gamma$. We compare the results with the outcome of extensive simulations performed with the reduced-rejection-rate Monte Carlo method Baldassi (2017), in which $\Gamma$ is initialized at $2.5$ and gradually brought down to [math] in regular small steps, at constant temperature, and fixing the total simulation time to $\tau Ny\cdot 10^{4}$ (as to keep constant the number of Monte Carlo sweeps when varying $N$ and $y$). The details are reported in the Appendix Sec. C. The size of the systems, the number of samples and the number of Trotter replicas are scaled up to large values so that both finite size effects and the quantum limit are kept under control. A key point is to observe that the results do not degrade with the number of Trotter replicas: the average ground state energy approaches a limiting value, close to the theoretical prediction, in the large $y$ quantum limit. The results appear to be rather insensitive to both $N$ and the simulation time scaling parameter $\tau$. This indicates that Monte Carlo appears to be able to equilibrate efficiently, in a constant (or almost constant) number of sweeps, at each $\Gamma$. The analytical prediction for the classical energy only appears to display a relatively small systematic offset (due to the static approximation) at intermediate values of $\Gamma$, while it is very precise at both large and small $\Gamma$; the expectation of the total Hamiltonian on the other hand is in excellent agreement with the simulations (see Appendix Sec. C).

In the same plot we display the behavior of classical SA simulated with a standard Metropolis-Hastings scheme, under an annealing protocol in $\beta$ that would follow the same theoretical curve as SQA if the system were able to equilibrate (see Appendix Sec. C): as expected Horner (1992), SA gets trapped at very high energies (increasing with problem size; in the thermodynamic limit it is expected that SA would remain stuck at the initial value $0.5N$ of the energy for times which scale exponentially with $N$). Alternative annealing protocols yield analogous results; the exponential scaling with $N$ of SA on binary perceptron models had also been observed experimentally in previous results, e.g. in refs. Baldassi et al. (2016b, 2007).

In the inset of Fig. 2 we report the analytical prediction for the transverse overlap parameter $q_{1}$, which quite remarkably reproduces fairly well the average overlap as measured from simulations.

In Fig. 3 we provide the profiles of the the classical energy minima found for different values of $\Gamma$ in the case of SQA and different temperatures for SA. These results are computed analytically by the cavity method (see Materials and Methods and SI for details) by evaluating which is the most probable energy found at a normalized Hamming distance $d$ from a given configuration. As it turns out, throughout the annealing process, SQA follows a path corresponding to wide valleys while SA gets stuck in steep metastable states. The quantum fluctuations reproduced by the SQA process drive the system to converge toward wide flat regions, in spite of the fact that they are exponentially rare compared to the narrow minima.

The physical interpretation of these results is that quantum fluctuations lower the energy of a cluster proportionally to its size or, in other words, that quantum fluctuations allow the system to lower its kinetic energy by delocalizing, see Refs. Foini et al. (2010); Markland et al. ; Biroli and Zamponi (2012) for related results. Along the process of reduction of the transverse field we do not observe any phase transition which could induce a critical slowing down of the quantum annealing process and we expect SQA and QA to behave similarly Bapst and Semerjian (2013); Bapst et al. (2013).

This is in agreement with the results of a direct comparison between the real time quantum dynamics and the SQA on small systems ($N=21$): as reported in the Appendix Sec. E, we have performed extensive numerical studies of properly selected small instances of the binary perceptron problem, comparing the results of SQA and QA and analyzing the results of the QA process and the properties of the Hamiltonian. To reproduce the conditions that are known to exist at large values of $N$, we have selected instances for which a fast annealing schedule SA gets trapped at some positive fraction of violated constraints, and yet the problems display a sufficiently high number of solutions. We found that the agreement between SQA and QA on each sample is excellent. The measurements on the final configurations reached by QA qualitatively confirm the scenario described above, that QA is attracted towards dense low-energy regions without getting stuck during the annealing process. Finally, the analysis of the gap between the ground state of the system and the first excited state as $\Gamma$ decreases shows no signs of the kind of phenomena which would typically hamper the performance of QA in other models: there are no vanishingly small gaps at finite $\Gamma$ (cf. the discussion in the introduction). We benchmarked all these results with “randomized” versions of the same samples, in which we randomly permuted the classical energies associated to each spin configuration, so as to keep the distribution of the classical energy levels while destroying the geometric structure of the states. Indeed, for these randomized samples, we found that the gaps nearly close at finite $\Gamma\simeq 0.4$, and that correspondingly the QA process fails to track the ground state of the system, resulting in a much reduced probability of finding a solution to the problem.

To reproduce the conditions that are known to exist at large values of $N$ we have selected instances for which a fast annealing schedule SA gets trapped at some positive fraction of violated constraints and yet the problems display a sufficiently high number of solutions. We have then compared the behavior of SQA and the real time quantum dynamics studied by the Lanczos method as discussed in Schneider et al. (2016). The agreement between SQA and QA is … almost perfect.

As concluding remarks we report that the models with $n=0$ and $n=1$ have phase diagrams which are qualitatively very similar (for the sake of simplicity, here we reported the $n=0$ case only). The former presents at very small positive values of $\Gamma$ a collapse of the density matrix onto the classical one whereas the latter ends up in the classical state only at $\Gamma=0$.

For the sake of completeness, we have checked that the performance of SQA in the $y\to\infty$ quantum limit extends to more complex architectures which include hidden layers; the details are reported in the Appendix Sec. D.2.

V Conclusions

We conclude by noticing that, at variance with other studies on spin glass models in which the evidence for QA outperforming classical annealing was limited to finite values of $y$, thereby just defining a different type of classical SA algorithms, in our case the quantum limit coincides with the optimal behavior of the algorithm itself. We believe that these results could play a role in many optimization problems in which optimality of the cost function needs to also meet robustness conditions (i.e. wide minima). As far as learning problems are concerned, it is worth mentioning that for the best performing artificial neural networks, the so called deep networks LeCun et al. (2015), there is numerical evidence for the existence of rare flat minima Keskar et al. (2016), and that all the effective algorithms always include effects of systematic injected noise in the learning phase Bottou et al. (2016), which implies that the equilibrium Gibbs measure is not the stationary measure of the learning protocols. For the sake of clarity we should remark that our results are aimed to suggest that QA can equilibrate efficiently whereas SA cannot, i.e. our notion of quantum speed up is relative to the same algorithmic scheme that runs on classical hardware. Other classical algorithms for the same class of problems, besides the above-mentioned ones based on the RE and the SQA itself, have been discovered Braunstein and Zecchina (2006); Baldassi et al. (2007); Baldassi (2009); Baldassi and Braunstein (2015); Hubara et al. (2016); however, all of these algorithms are qualitatively different from QA, which can provide a huge speed up by manipulating single bits in parallel. Thus, the overall solving time in a physical QA implementation (neglecting any other technological considerations) would have, at worst, only a mild dependence on $N$.

Our results provide further evidence that learning can be achieved through different types of correlated fluctuations, among which quantum tunneling could be a relevant example for physical devices.

Acknowledgements.

The authors thank G. Santoro, B. Kappen and F. Becca for discussions.

Appendix A Theoretical analysis by the replica method

We present here the analytical calculations performed to derive all the theoretical results mentioned in the main text. For completeness, we report all the relevant formulas and definitions here, even those that were already introduced in the main text.

The Hamiltonian operator of a model of $N$ quantum spins with an energy term acting in the longitudinal direction $z$ and a magnetic field $\Gamma$ acting in the transverse direction $x$ is written as:

[TABLE]

where $\hat{\sigma}_{j}^{z}$ and $\hat{\sigma_{j}}^{x}$ are the spin operators (Pauli matrices) in the $z$ and $x$ directions. We want to study the partition function:

[TABLE]

By using the Suzuki-Trotter transformation, we end up with a classical effective Hamiltonian acting on a system of $y$ interacting Trotter replicas, to be studied in the limit $y\to\infty$:

[TABLE]

where the $\sigma_{j}^{a}=\pm 1$ are Ising spins, $a\in\left\{1,\dots,y\right\}$ is a replica index with periodic boundary conditions $\sigma_{j}^{y+1}\equiv\sigma_{j}^{1}$, and we have defined:

[TABLE]

In the following, we will just use $\sigma^{a}$ to denote the configuration of one Trotter replica, $\left\{\sigma_{j}^{a}\right\}_{j}$; we will always use the indices $a$ or $b$ for the Trotter replicas and assume that they range in $1,\dots,y$; we will also use $j$ for the site index and assume that it ranges in $1,\dots,N$.

The effective partition function for a given $y$ reads:

[TABLE]

Here, we first study the binary perceptron case in which the longitudinal energy $E$ is defined in terms of a set of $\alpha N$ patterns $\left\{\xi^{\mu}\right\}_{\mu}$ with $\mu\in\left\{1,\dots,\alpha N\right\}$, where each pattern is a binary vector of length $N$, $\xi_{j}^{\mu}=\pm 1$:

[TABLE]

where $\Theta\left(x\right)$ is the Heaviside step function: $\Theta\left(x\right)=1$ if $x>0$, $\Theta\left(x\right)=0$ otherwise. The energy thus simply counts the number of classification errors of the perceptron, assuming that the desired output for each pattern in the set is $1$ (this choice can always be made without loss of generality for this model when the input patterns are random i.i.d. as described below). A different form for the energy function is treated in sec. A.1.

We consider the case in which the patterns entries are extracted randomly and independently from an unbiased distribution, $P\left(\xi_{j}^{\mu}\right)=\frac{1}{2}\delta\left(\xi_{j}^{\mu}-1\right)+\frac{1}{2}\delta\left(\xi_{j}^{\mu}+1\right)$, and we want to study the typical properties of this system by averaging over the quenched disorder introduced by the patterns. We use the replica method, which exploits the transformation:

[TABLE]

where $\left\langle\cdot\right\rangle_{\xi}$ denotes the average over the disorder. We thus need to replicate the whole system $n$ times, and therefore we have two replica indices for each spin. We will use indices $c,d=1,\dots,n$ for the “virtual” replicas introduced by the replica method,111Note that the parameter $n$ has a different meaning in main text, cf. sec. A.1. to distinguish them from the indices $a$ and $b$ used for the Trotter replicas. The average replicated partition function of eq. (10) is thus written as:

[TABLE]

where we changed the sum over all configurations into an ($n\times y\times N$-dimensional) integral, using the customary notation $d\mu\left(\sigma\right)=\delta\left(\sigma-1\right)+\delta\left(\sigma+1\right)$ with $\delta\left(\cdot\right)$ denoting the Dirac-delta distribution. Here and in the following, all integrals will be assumed to range over the whole $\mathbb{R}$ unless otherwise specified.

We introduce new auxiliary variables $\lambda_{\mu}^{ca}=\frac{1}{\sqrt{N}}\sum_{j}\xi_{j}^{\mu}\sigma_{j}^{ca}$ via additional Dirac-deltas:

[TABLE]

We then use the integral representation of the delta $\delta\left(x\right)=\int\frac{d\hat{x}}{2\pi}e^{ix\hat{x}}$, and perform the average over the disorder, to the leading order in $N$:

[TABLE]

Next, we introduce the overlaps $q^{ca,db}=\frac{1}{N}\sum_{j}\sigma_{j}^{ca}\sigma_{j}^{db}$ via Dirac-deltas (note that due to symmetries and the fact that the self-overlaps are always $1$ we have $ny\left(ny-1\right)/2$ overlaps overall), expand those deltas introducing conjugate parameters $\hat{q}^{ca,db}$ (as usual for these parameters in these models, we absorb away a factor $i$ and integrate them along the imaginary axis, without explicitly noting this), and finally factorize over the site and pattern indices:

[TABLE]

We now introduce the replica-symmetric (RS) ansatz for the overlaps:

[TABLE]

and analogous for the conjugate parameters $\hat{q}^{ca,db}$.

Note that this is the so-called “static approximation” since we neglect the dependency of the overlap from the distance along the Trotter dimension; however, we have kept the interaction term $\gamma\sum_{ca}\sigma^{ca}\sigma^{c\left(a+1\right)}$ and inserted it in the $G_{S}$ term (rather than writing it in terms of the overlap $q^{ca,c\left(a+1\right)}$ and inserting it in the $G_{E}$ term where it would have been rewritten as $\gamma q_{1}$). This difference, despite its inconsistency, is the standard procedure when performing the static approximation, and is justified a posteriori from the comparison with the numerical simulation results. We obtain:

[TABLE]

The entropic term $G_{S}$ can be explicitly computed as

[TABLE]

where the notation $Dz=dz\,\frac{1}{\sqrt{2\pi}}e^{-\frac{x^{2}}{2}}$ is a shorthand for a Gaussian integral, and we used twice the Hubbard-Stratonovich transformation $e^{\frac{1}{2}b}=\int Dz\,e^{z\sqrt{b}}$. The expression between square brackets in the last line is the partition function of a 1-dimensional Ising model of size $y$ with uniform interactions $J=\gamma$ and uniform fields $h=z_{1}\sqrt{\hat{q}_{1}-\hat{q}_{0}}+z_{0}\sqrt{\hat{q}_{0}}$ and can be computed by the well-known transfer matrix method. Note however that while usually in the analysis of the 1D Ising spin model it is sufficient to keep the largest eigenvalue of the transfer matrix in the thermodynamic limit $y\to\infty$, in this case instead we need to keep both eigenvalues, since the interaction term scales with the size of the system. The result is:

[TABLE]

In the limit of small $n$ we obtain:

[TABLE]

Note that in the limit of large $y$ the term $\gamma y$ tends to$-K$ up to terms of order $y^{-1}$.

The energetic term $G_{E}$ is computed similarly, by first performing two Hubbard-Stratonovich transformations which allow to factorize the indices $c$ and $a$, and then explicitly performing the inner integrals:

[TABLE]

where $H\left(x\right)=\frac{1}{2}\mathrm{erfc}\left(\frac{x}{\sqrt{2}}\right)$. In the limit of small $n$ and of large $y$ we finally obtain:

[TABLE]

Using equations (27) and (29), we obtain the expression for the action:

[TABLE]

In order to obtain a finite result in the limit of $y\to\infty$, we assume the following scalings for the conjugated order parameters:

[TABLE]

With these, we find the following final expressions:

[TABLE]

The parameters $q_{0}$, $q_{1}$, $\hat{p}_{0}$ and $\hat{p}_{1}$ are found by solving the system of equations obtained by setting the partial derivatives of $\phi$ with respect to those parameters to [math]:

[TABLE]

Once these are found, we can use them to compute the action $\phi$ and the average values of the longitudinal energy and the transverse fields, and finally of the Hamiltonian:

[TABLE]

where the notation $\overline{\left\langle\cdot\right\rangle_{\xi}}$ denotes the fact that we performed both the average over the quenched disorder and the thermal average.

A.0.1 Small $\Gamma$ limit

It can be verified that in the limit $\Gamma\to 0$ the equations (33)-(37) reduce to the classical case, in the RS description. In this limit, $q_{1}\to 1$ (i.e., the Trotter replicas collapse), which leads to:

[TABLE]

For $\Gamma=0$ and $q_{1}=1$ we also have the identity:222This follows from $\int Dz_{1}\cosh\left(a\,z_{1}+b\,z_{0}\right)=e^{\frac{a^{2}}{2}}\cosh\left(b\,z_{0}\right)$.

[TABLE]

Putting these two expressions back in eq. (33) we recover the classical expression where $\hat{p}_{0}$ assumes the role of the usual conjugate parameter $\hat{q}$ in the RS analysis of ref. Krauth and Mézard (1989).

In order to study in detail how this classical limit is reached, however, we need to expand the saddle point equations around this limit. To to this, we define $\epsilon=1-q_{1}\ll 1$. From equation (A), expanding to the leading order, we obtain the scaling $\hat{p}_{1}=\frac{\hat{c}_{1}}{\sqrt{\epsilon}}$, with

[TABLE]

Then, we use this scaling in equation (A) and we expand it, first using $\beta\Gamma\ll 1$ and then $\epsilon\ll 1$. We obtain the approximate expression:

[TABLE]

where $F\left(x\right)=\frac{\sqrt{\pi}}{2}e^{-x^{2}}\mathrm{erfi}\left(x\right)$ is the Dawson’s function. For a given $\beta$ (from which we obtain $\hat{c}_{1}$ via eq. (47)), this equation can be solved numerically to obtain $\epsilon$ (and thus $q_{1}$ and $\hat{p}_{1}$) as a function of $\Gamma$. This expression has always the solution $\epsilon=0$, which correspond to the purely classical case. There is a critical $\Gamma$ below which $\epsilon=0$ is also the only solution; above that, two additional solutions appear at $\epsilon>0$, of which the largest is the physical one. Therefore, the classical limit is not achieved continuously, but rather with a first-order transition (although the step is tiny).

A.1 Energy function with stability

We can generalize the energy function eq. (11) to take into account, for those patterns that are misclassified, by how much the classification is wrong:

[TABLE]

The previous case is recovered by setting $r=0$. Here, we study the case $r=1$. Note that this parameter is called $n$ in the main text: that notation was borrowed from ref. Horner (1992), but here we change it in order to avoid confusion with the number of replicas. While the ground states in the SAT phase of the classical model are unaffected, the system can have different properties for finite $\beta$.

This change only affects the $\mathcal{G}_{E}$ term. Equation (28) becomes (with the definition of eq. (37)):

[TABLE]

In the limit of large $y$ we have the modified version of eq. (29):

[TABLE]

The saddle point equations (38) and (A) become:

[TABLE]

where

[TABLE]

A.1.1 Small $\Gamma$ limit

As in the previous case, it can be checked that for $\Gamma\to 0$, we have $q_{1}\to 1$ and eq. (51) becomes the expression for the classical model under the RS ansatz:

[TABLE]

where $k_{0}\left(z_{0}\right)=z_{0}\sqrt{\frac{q_{0}}{1-q_{0}}}$. Also, eq. (46) still holds, and $\hat{p}_{0}$ takes the role of the usual parameter $\hat{q}$ in the classical RS analysis. In this case, however, we no longer have $\hat{p}_{1}\to\infty$; rather, it tends to a finite value:

[TABLE]

Therefore, the scaling of $\epsilon=1-q_{1}$ is different in this case. We find (using the definition of eq. (35)):

[TABLE]

Therefore, the convergence to the classical case is smooth.

Appendix B Estimation of the local energy and entropy landscapes with the cavity

method

In order to compute the local landscapes of the energy and the entropy around a reference configuration (Fig. 3), we used the Belief Propagation (BP) algorithm, a cavity method message-passing algorithm that has been successfully employed numerous times for the study of disordered systems Mézard and Montanari (2009). In the case of single-layer binary perceptrons trained on random unbiased i.i.d. patterns, it is believed that the results of this algorithm are exact in the limit of $N\to\infty$, at least up to the critical value $\alpha_{c}\approx 0.83$ Mézard (1989).

For a full explanation of the BP equations for binary perceptrons, we refer the interested reader to the Appendix of ref. Baldassi et al. (2016b). Here, we provide only a summary. The BP equations involve two sets of quantities (called “messages”), representing cavity marginal probabilities associated with each edge in a factor graph representation of the (classical) Boltzmann distribution induced by the energy function (11). To each edge in the graph linking the variable node $i$ with the factor node $\mu$, are associated two messages, $m_{i\to\mu}$ and $\hat{m}_{\mu\to i}$. These are determined by solving iteratively the following system of equations:

[TABLE]

where:

[TABLE]

(as for the previous section, we used the definition $H\left(x\right)=\frac{1}{2}\mathrm{erfc}\left(\frac{x}{\sqrt{2}}\right)$.)

Once a self-consistent solution is found, these quantities can be used to compute, using standard formulas, all thermodynamic quantities of interest, in particular the typical (equilibrium) energy and the entropy of the system. A numerically accurate implementation of these equations is available at ref. FBP .

It is also possible to compute those same thermodynamic quantities in a neighborhood of some arbitrary reference configuration $w=\left\{w_{i}\right\}_{i}$. This is achieved by adding an external field in the direction of that configuration, which amounts at this simple modification of eq. (57):

[TABLE]

By varying the auxiliary parameter $\lambda$, we can control the size of the neighborhood under consideration (the larger $\lambda$, the narrower the neighborhood); the typical normalized Hamming distance from the reference of the configurations that are considered by this modified measure can be obtained from the fixed-point BP messages for any given $\lambda$ by this formula:

[TABLE]

where the $m_{i}$ are the total magnetizations:

[TABLE]

In order to produce the energy landscape plots of Figs. 3a and 3b, we simply ran this algorithm at infinite temperature, varying $\lambda$ and plotting the energy density shift from the center as a function of $d$. This gives us an estimate of the most probable energy density shift which would be obtained by moving in a random point at distance $d$ from the reference.

The plot in Fig. 3c was similarly obtained by setting the temperature to [math] and computing the entropy density instead, which in this context is then simply the natural logarithm of the number of solutions in the given neighborhood, divided by $N$.

Appendix C Numerical simulations details of the annealing protocols

C.1 Quantum annealing protocol

In this section we provide the details of the QA results presented in Fig. 2. The simulations were performed using the RRR Monte Carlo method Baldassi (2017). We fixed the total number of spin flip attempts at $\tau Ny\cdot 10^{4}$ and followed a linear protocol for the annealing of $\Gamma$, starting from $\Gamma_{0}=2.5$ and reaching down $\Gamma_{1}=0$. We actually divided the annealing in $30\tau$ steps, where during each step $\Gamma$ was kept constant and decreased by $\Delta\Gamma=\frac{\Gamma_{0}-\Gamma_{1}}{30\tau}$ after each step. In the figure, we have shown the results for $N=4001$ and $\tau=4$; the results for $N=1001,2001$ and for $\tau=1,2$ were essentially indistinguishable at that level of detail.

C.2 Classical simulated annealing protocol

The results for SA presented in Fig. 2 used an annealing protocol in $\beta$ designed to make a direct comparison to QA: we found analytically a curve $\beta_{\mathrm{equiv}}\left(\Gamma\right)$ such that the classical equilibrium energy would be equal to the longitudinal component of the quantum system energy, eq. (43). The classical equilibrium energy was computed from the equations in ref. Krauth and Mézard (1989). The result is shown in Fig. 4. The vertical jump to $\beta=20$ is due to the transition mentioned in sec. A.0.1; as shown in Fig. 2, the SA protocol in the regime we tested gets stuck well before this transition.

The SA annealing protocol thus consisted in setting $\beta=\beta_{\mathrm{equiv}}\left(\Gamma\right)$ and decreasing linearly $\Gamma$ from $2.5$ to [math], like for the QA case. We fixed the total number of spin flip attempts at $\tau N\cdot 10^{4}$ and used $\tau=4,8,16$; as for the QA case, the annealing process was divided in $30\tau$ steps.

Other more standard annealing protocols (e.g. linear or exponential or logarithmic) yielded very similar qualitative results, as expected from the analysis of ref. Horner (1992).

Appendix D Additional numerical results on the annealing processes

D.1 Additional comparisons between theory and simulations

Fig. 2 compares the result of Monte Carlo simulations with the theoretical predictions for the classical component of the energy, eq. (43), and the transverse overlap, eq. (A). Fig. 5 compares the same simulation results with the analytical curves at finite $y$ instead. This shows a relatively small systematic offset (due to the static approximation) at intermediate values of $\Gamma$, while the agreement is good at both large and small $\Gamma$.

Fig. 6 shows the comparison with the $y\to\infty$ curve for the expectation of the full quantum Hamiltonian, eq. (42), using the same data. The agreement is remarkable, and a close inspection reveals that the curves from the simulation tend towards the theoretical one as $y$ increases, i.e. in the quantum limit.

D.2 Experiments with two-layer networks

We performed additional experiments using two-layer fully-connected binary networks, the so-called committee machines. Previous results obtained with the robust-ensemble measure Baldassi et al. (2016a) showed that this case is quite similar to that of single layer networks. In particular, standard Simulated Annealing suffers from an exponential slow-down as the system size increases even moderately, while algorithms that are able to target the dense states do not suffer from the trapping in meta-stable states. Indeed, we found the latter feature to be true in the quantum annealing scenario.

The model in this case is defined by a modified energy function (cf. eq. (11)):

[TABLE]

where now the $N$ spin variables are divided in groups of $K$ hidden units, and consequently the spin variables $\sigma_{kj}$ have two indices, $k=1,\dots,K$ for the hidden unit and $j=1,\dots,N/K$ for the input. Notice that the input size is reduced $K$-fold with respect to the previous case. The output of these machines is simply decided by the majority of the outputs of the individual units, and the energy still counts the number of errors. The Suzuki-Trotter transformation proceeds in exactly the same way as for the previous cases.

Like for the single-layer case, we tested the case of $\alpha=0.4$ at $\beta=20$, and we used $K=5$ units. We tested different values of $N=1005,2005,4005$ with different values of the Trotter replicas $y=32,64,128$ (only $y=32$ for $N=4005$) at a fixed overall running time of $yN\tau\cdot 10^{4}$ spin flip attempts, with $\tau=4$ (cf. Fig. 2). The MC algorithm and the annealing protocols were also unchanged. The results are shown in Fig. 7: all these tests produce curves which are almost indistinguishable at this level of detail for different $N$, and that seemingly tend to converge to some limit curve for increasing $y$ (while being almost overlapping at small transverse field $\Gamma$), consistently with the single-layer scenario.

Appendix E Real-time Quantum Annealing on small samples

E.1 Numerical methods

Computing the evolution of the system under Quantum Annealing amounts at solving the time-dependent Schr dinger equation for the system

[TABLE]

where we set $\hslash=1$ for simplicity. In our case, the time dependence of the Hamiltonian $H$ comes in through the varying transverse magnetic field $\Gamma\left(t\right)$. We assume that $\Gamma$ varies linearly with time between some starting value $\Gamma_{0}$ and [math], in a total time $t_{\max}$. Therefore, the final Hamiltonian is reduced to the purely classical case, $\hat{H}\left(t_{\max}\right)=E$.

In the following, we will always work in the basis of the final Hamiltonian, in which every eigenvector $\left|\sigma\right\rangle$ corresponds to a configuration $\sigma\in\left\{-1,+1\right\}^{N}$ of the spins in the $z$ direction. Therefore, we represent $\left|\psi\left(t\right)\right\rangle$ with a complex-valued vector of length $2^{N}$ with entries $\left\langle\sigma|\psi\left(t\right)\right\rangle$; similarly, the $\hat{H}\left(t\right)$ operator is represented by a matrix of size $2^{N}\times 2^{N}$, $H\left(\sigma,\sigma^{\prime}\right)=\left\langle\sigma\left|\hat{H}\right|\sigma^{\prime}\right\rangle$. The structure of this matrix is very sparse: the diagonal elements $H\left(\sigma,\sigma\right)$ correspond to the classical energies $E\left(\sigma\right)$, while the only non-zero diagonal elements are those elements $H\left(\sigma,\sigma^{\prime}\right)$ such that $\sigma$ and $\sigma^{\prime}$ are related by a single spin flip, in which case the value is $-\Gamma$.

In our simulations, the initial state $\left|\psi\left(0\right)\right\rangle$ was set to the ground state of the system at $\Gamma\to\infty$, i.e. with all the spins aligned in the $x$ direction; in our basis, this corresponds to a uniform vector, $\left\langle\sigma|\psi\left(0\right)\right\rangle=N^{-\nicefrac{{1}}{{2}}}$ for all $\sigma$. We simulated the evolution of the system by the short iterative Lanczos (SIL) method Schneider et al. (2016): we compute the evolution at fixed $\Gamma$ for a short time interval $\Delta t$, then lower $\Gamma$ by a small fixed amount $\Delta\Gamma$, and iterate. The total evolution time is thus $t_{\max}=\frac{\Gamma_{0}}{\Delta\Gamma}\Delta t$. Numerical accuracy can be verified by scaling both these steps by a fixed amount and observing no significant difference in the outcome. The evolution is computed by the Lanczos algorithm with enough iterative steps to ensure sufficient accuracy, as determined by observing that increasing the number of steps does not change the outcomes significantly. In the simulations presented here, we set $\Gamma_{0}=5$, $\Delta\Gamma=10^{-3}$ and $\Delta t=0.2$, and we used $10$ steps in the Lanczos iterations.

At the end of the annealing process, we could retrieve the final probability distribution for each configuration of the spins as $p\left(\sigma\right)=\left|\left\langle\sigma|\psi\left(t_{\max}\right)\right\rangle\right|^{2}$.

E.2 Sample selection

Given the exponential scaling with $N$ of the SIL algorithm, simulations are necessarily restricted to small values of $N$. We used $N=21$. At these system sizes, there is a very large sample-to-sample variability. Furthermore, the energy barriers are generally small enough for the classical Simulated Annealing to perform well.

In order to obtain small but challenging samples, in which we could also study the structure of the solutions, we proceeded as follows: we extracted at random $450$ samples with $P=17$ patterns each (corresponding to $\alpha\simeq 0.81$, close to the critical value of $0.83$ which is valid for large systems), and selected those which had at least a certain minimum number of solutions (note that in such small systems we can easily enumerate all of the $2^{21}\simeq 2\cdot 10^{6}$ configurations and check their energy). We arbitrarily chose $21$ solutions as the threshold. We then ran both Simulated Annealing with a fast schedule (with $\tau=1)$ and Simulated Quantum Annealing with $\tau=1$ and a large number of Trotter replicas ($y=512$), and selected those samples in which SA failed while SQA succeeded. This left us with $20$ samples, which we then analyzed in detail and over which we performed the real-time QA simulations.

E.2.1 Randomized samples

For each of the selected samples, we generated a corresponding randomized version by permuting randomly the values of the energy associated to each configuration. This procedure maintains unaltered the spectrum of the energies (so that for example the classical Boltzmann distribution at thermodynamic equilibrium remains unchanged), but completely destroys the geometric features of the energy landscape. We used these randomized samples as a benchmark against the measurements performed in our analysis.

E.3 Analysis

E.3.1 QA vs SQA

We compared the results of real-time QA with the SQA Monte Carlo results, analyzing each of the $20$ selected samples individually. In particular, we compared the values of the average longitudinal energy as a function of $\Gamma$ for the two algorithms. As shown in Fig. 8, the agreement is excellent, and the system always gets very close to zero energy. In the same figure, we show that the same annealing protocol however gives substantially different (and rather worse) results on the randomized samples, reflecting the fact that the geometrical features of the landscape are crucial (we verified on a few cases that the results on the randomized samples could be improved by slowing down the annealing process, but we could not get to the same results as for the original systems even with a $100$-fold increase in total time). Note that the sample-to-sample variability in these curves appears to be fairly small due to our sample filtering process; we verified in a preliminary analysis that the agreement is generally excellent also without the filtering conditions, e.g. on instances that have no solutions at all.

E.3.2 Other measurements

We also performed a number of measurements on the final configuration reached by the QA algorithm (both for the original samples and the randomized ones) and studied the properties of the final probability distribution $p\left(\sigma\right)$. These are the quantities that we computed, reported in table 1:

• •

The average value of the energy $\left\langle E\right\rangle=\sum_{\sigma}E\left(\sigma\right)p\left(\sigma\right)$.

• •

The probability of finding a solution $P_{\mathrm{SOL}}=\sum_{\sigma:E\left(\sigma\right)=0}p\left(\sigma\right)$.

• •

The probability and the energy of the most probable configuration, $p\left(\sigma^{\star}\right)$ and $E\left(\sigma^{\star}\right)$, where $\sigma^{\star}=\arg\max_{\sigma}p\left(\sigma^{\star}\right)$.

• •

The inverse participation ratio $\mathrm{IPR}=\sum_{\sigma}p\left(\sigma\right)^{2}$, to assess the concentration of the final distribution. (Qualitatively analogous results are obtained using the Shannon entropy.) This measure however does not take into account the geometric structure of the distribution: for instance, if $p\left(\sigma\right)$ were non-zero on just to configurations $\sigma^{1}$ and $\sigma^{2}$, the $\mathrm{IPR}$ would be very high, but it would not be able to discriminate between the cases in which $\sigma^{1}$and $\sigma^{2}$ are close to each other or far apart in Hamming distance.

• •

The mean distance between configurations, defined as $\bar{d}=\sum_{\sigma,\sigma^{\prime}}p\left(\sigma\right)p\left(\sigma^{\prime}\right)d\left(\sigma,\sigma^{\prime}\right)$, where $d\left(\sigma,\sigma^{\prime}\right)$ is the normalized Hamming distance between configurations. This measure is useful since it reflects the geometric features of the final measure: it can only be low if the mass of the probability is concentrated spatially (in particular, it is zero if and only if $p\left(\sigma\right)$ is a delta function).

As can be seen from the table, the results are generally in agreement with the qualitative picture described in the main text, especially when compared to the randomized benchmark: the system is able to reach very low energies $\left\langle E\right\rangle$, the probability of solving the problem $P_{\mathrm{SOL}}$ is very high, the measure is rather concentrated on a few good configurations and those configurations are close to each other (high $\mathrm{IPR}$, low $\bar{d}$).

For the original samples only, we also looked at the final configuration from the Monte Carlo SQA process, $\sigma_{\mathrm{SQA}}$, and computed its ranking according to $p\left(\sigma\right)$, which we denoted as $r_{\mathrm{SQA}}$. A ranking of $1$ implies $\sigma_{\mathrm{SQA}}=\sigma^{\star}$. All rankings are very small, the largest ones generally corresponding to samples with the largest number of solutions and the less concentrated distributions. This further attests to the good agreement between the QA and SQA processes.

E.3.3 Local entropies

In order to assess whether the denser ground states were favored in the final configuration with respect to more isolated solutions we compared the mean local entropy curves weighted according to $p\left(\sigma\right)$ with those averaged over all the solutions. More precisely, we define $C\left(n\right)$ as the set of the $n$ configurations with highest probability, and $n_{w}$ as the number of configurations required to achieve a cumulative probability of $t$, i.e. the lowest $n$ such that $\sum_{\sigma\in C\left(n\right)}p\left(\sigma\right)\geq w$. We also define $K\left(\sigma,d\right)$ as the number of solutions at normalized Hamming distance from $\sigma$ lower or equal to $d$. Then the mean local entropy curve weighted with $p$ is then defined as:

[TABLE]

Denoting by $\mathcal{S}=\left\{\sigma|E\left(\sigma\right)=0\right\}$ the set of all the solutions, we also compute the flat average of the local entropies over $\mathcal{S}$:

[TABLE]

If $p$ concentrates on denser solutions, we expect that the $\phi_{w}$ curves should be generally higher than the $\phi_{\mathrm{SOL}}$ curves. Indeed, the results confirm this scenario, as shown in Fig. 9, where we used $w=0.9$. (This value ensured that $C\left(n_{w}\right)\subseteq\mathcal{S}$ for all samples and thus that all the local entropies are finite; apart from this, the results are quite insensitive to the choice of $w$.) Note that, in the limit of large system sizes, the $\phi_{\mathrm{SOL}}$ curves would be dominated by isolated solutions and display a gap around zero distances; the fact that this is not visible in Fig. 9 is purely a finite size effect; the $\phi_{w}$ curves on the other hand should be roughly comparable to those shown in Fig. 3.

E.3.4 Energy gaps

As mentioned in the introduction of the main text, it is well known that, according to the adiabatic theorem, the effectiveness of the QA process depends on the relation between the rate of change of the Hamiltonian and the size of the gap between the ground state of the system $H_{0}$ and the first excited state $H_{1}$: smaller gaps require a slower annealing process. Therefore, we performed a static analysis of the energy spectrum of each of the $20$ samples at varying $\Gamma$, and computed the gap $H_{1}-H_{0}$, comparing the results with those for the randomized versions of the samples. The results are shown in Fig. 10. For the original samples, the gap only vanishes in the limit of $\Gamma\to 0$ (which is expected since the ground state at $\Gamma=0$ is degenerate). For the randomized samples, on the other hand, the gap nearly closes at non-zero $\Gamma$, displaying the characteristics of an “avoided crossing” (see the figure upper inset), which is the type of phenomenon that is known to hamper the performance of QA algorithms. Indeed, the values of $\Gamma$ where these avoided crossings occur are precisely those at which the mean value of $H$ found by the QA algorithm deviates from the ground state $H_{0}$, thereby getting stuck as shown in Fig. 8.