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 scienceMoore 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 stateBorn 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 glassesAltshuler 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 wholeMacKay (2003)
. In this respect, it is worth mentioning that in deep learningLeCun 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 problemKrauth and Mézard (1989)
. This network maps vectors ofinputs to binary outputs through the non linear function , where is the vector of synaptic weights. Given input patterns with and their corresponding desired outputs , the learning problem consists in finding
such that all input patterns are simultaneously classified correctly, i.e.for all . Both the components of the input vectors and the outputs
are independent identically distributed unbiased random variables (). 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 of the binary perceptron is proportional to the number of classification errors and can be written as
where is the Heaviside step function: if , otherwise. When the argument of the function is positive, the perceptron is implementing the wrong input-output mapping. The exponent 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 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 the problem is computationally hard for SA processes Horner (1992): in the large limit, a detailed balanced stochastic search process gets stuck in metastable states at energy levels of order above the ground states.
Following the standard SQA approach, we identify the binary variableswith one of the components of physical quantum spins, say , and we introduce the Hamiltonian operator of a model of quantum spins with the perceptron term of Eq. (1) acting in the longitudinal direction and a magnetic field acting in the transverse direction . The interpolating Hamiltonian reads:
where and are the spin operators (Pauli matrices) in the and directions. For one recovers the classical optimization problem. The QA procedure consists in initializing the system at large and , and slowly decreasing to . To analyze the low temperature phase diagram of the model we need to study the average of the logarithm of the partition function This can be done using the Suzuki-Trotter transformation which leads to the study of a classical effective Hamiltonian acting on a system of interacting Trotter replicas of the original classical system coupled in an extra dimension:
where the are Ising spins, is a replica index with periodic boundary conditions , and
The replicated system needs to be studied in the limit 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 diverges as the transverse field goes to . 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, i.e. a “local free entropy” where is a Lagrange parameter that controls the extensive size of the region around a reference configuration . One can then build a new Gibbs distribution , where has the role of an energy and of an inverse temperature: in the limit of large , this distribution concentrates on the maxima of . Upon restricting the values of to be integer (and large), takes a factorized form yielding a replicated probability measure where the effective energy is given by
As in the Suzuki-Trotter formalism, corresponds to a system with an overall energy given by the sum of individual “real replica energies” plus a geometric coupling term; in this case however the replicas interact with the “reference” configurations 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, and 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 , 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 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 to represent the overlaps along the Trotter dimension, . Although this approximation crudely neglects the dependency of from , 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 . We compare the results with the outcome of extensive simulations performed with the reduced-rejection-rate Monte Carlo method Baldassi (2017), in which is initialized at and gradually brought down to in regular small steps, at constant temperature, and fixing the total simulation time to (as to keep constant the number of Monte Carlo sweeps when varying and ). 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 quantum limit. The results appear to be rather insensitive to both and the simulation time scaling parameter . This indicates that Monte Carlo appears to be able to equilibrate efficiently, in a constant (or almost constant) number of sweeps, at each . 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 , while it is very precise at both large and small ; 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 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 of the energy for times which scale exponentially with ). Alternative annealing protocols yield analogous results; the exponential scaling with 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 , 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 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 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 (): 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 , 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 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 (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 , 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 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 and have phase diagrams which are qualitatively very similar (for the sake of simplicity, here we reported the case only). The former presents at very small positive values of a collapse of the density matrix onto the classical one whereas the latter ends up in the classical state only at .
For the sake of completeness, we have checked that the performance of SQA in the quantum limit extends to more complex architectures which include hidden layers; the details are reported in the Appendix Sec. D.2.
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, 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 .
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 quantum spins with an energy term acting in the longitudinal direction and a magnetic field acting in the transverse direction is written as:
where and are the spin operators (Pauli matrices) in the and directions. We want to study the partition function:
By using the Suzuki-Trotter transformation, we end up with a classical effective Hamiltonian acting on a system of interacting Trotter replicas, to be studied in the limit :
where the are Ising spins, is a replica index with periodic boundary conditions , and we have defined:
In the following, we will just use to denote the configuration of one Trotter replica, ; we will always use the indices or for the Trotter replicas and assume that they range in ; we will also use for the site index and assume that it ranges in .
The effective partition function for a given reads:
Here, we first study the binary perceptron case in which the longitudinal energy is defined in terms of a set of patterns with , where each pattern is a binary vector of length , :
where is the Heaviside step function: if , 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 (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, , 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:
where denotes the average over the disorder. We thus need to replicate the whole system times, and therefore we have two replica indices for each spin. We will use indices for the “virtual” replicas introduced by the replica method,111Note that the parameter has a different meaning in main text, cf. sec. A.1. to distinguish them from the indices and used for the Trotter replicas. The average replicated partition function of eq. (10) is thus written as:
where we changed the sum over all configurations into an (-dimensional) integral, using the customary notation with denoting the Dirac-delta distribution. Here and in the following, all integrals will be assumed to range over the whole unless otherwise specified.
We introduce new auxiliary variables via additional Dirac-deltas:
We then use the integral representation of the delta , and perform the average over the disorder, to the leading order in :
Next, we introduce the overlaps via Dirac-deltas (note that due to symmetries and the fact that the self-overlaps are always we have overlaps overall), expand those deltas introducing conjugate parameters (as usual for these parameters in these models, we absorb away a factor and integrate them along the imaginary axis, without explicitly noting this), and finally factorize over the site and pattern indices:
We now introduce the replica-symmetric (RS) ansatz for the overlaps:
and analogous for the conjugate parameters .
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 and inserted it in the term (rather than writing it in terms of the overlap and inserting it in the term where it would have been rewritten as ). 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:
The entropic term can be explicitly computed as
where the notation is a shorthand for a Gaussian integral, and we used twice the Hubbard-Stratonovich transformation . The expression between square brackets in the last line is the partition function of a 1-dimensional Ising model of size with uniform interactions and uniform fields
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, in this case instead we need to keep both eigenvalues, since the interaction term scales with the size of the system. The result is:
In the limit of small we obtain:
Note that in the limit of large the term tends to up to terms of order .
The energetic term is computed similarly, by first performing two Hubbard-Stratonovich transformations which allow to factorize the indices and , and then explicitly performing the inner integrals:
where . In the limit of small and of large we finally obtain:
In order to obtain a finite result in the limit of , we assume the following scalings for the conjugated order parameters:
With these, we find the following final expressions:
The parameters , , and are found by solving the system of equations obtained by setting the partial derivatives of with respect to those parameters to :
Once these are found, we can use them to compute the action and the average values of the longitudinal energy and the transverse fields, and finally of the Hamiltonian:
where the notation denotes the fact that we performed both the average over the quenched disorder and the thermal average.
a.0.1 Small limit
For and we also have the identity:222This follows from .
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 . From equation (A), expanding to the leading order, we obtain the scaling , with