## I Introduction

There is increasing evidence that artificial neural networks perform exceptionally well in complex recognition tasks lecun2015deep . In spite of huge numbers of parameters and strong non-linearities, learning often occurs without getting trapped in local minima with poor prediction performance ngiam2011optimization

. The remarkable output of these models has created unprecedented opportunities for machine learning in a host of applications. However, these practical successes have been guided by intuition and experiments, while obtaining a complete theoretical understanding of why these techniques work seems currently out of reach, due to the inherent complexity of the problem. In other words: in practical applications, large and complex architectures are trained on big and rich datasets using an array of heuristic improvements over basic stochastic gradient descent. These heuristic enhancements over a stochastic process have the general purpose of improving the convergence and robustness properties (and therefore the generalization properties) of the networks, with respect to what would be achieved with a pure gradient descent on a cost function.

There are many parallels between the studies of algorithmic stochastic processes and out-of-equilibrium processes in complex systems. Examples include jamming processes in physics, local search algorithms for optimization and inference problems in computer science, regulatory processes in biological and social sciences, and learning processes in real neural networks (see e.g. charbonneau2014fractal ; ricci2009cavity ; bressloff2014stochastic ; easley2010networks ; holtmaat2009experience

). In all these problems, the underlying stochastic dynamics are not guaranteed to reach states described by an equilibrium probability measure, as would occur for ergodic statistical physics systems. Sets of configurations which are quite atypical for certain classes of algorithmic processes become typical for other processes. While this fact is unsurprising in a general context, it is unexpected and potentially quite significant when sets of relevant configurations that are typically inaccessible for a broad class of search algorithms become extremely attractive to other algorithms.

Here we discuss how this phenomenon emerges in learning in large-scale
neural networks with low precision synaptic weights. We further show
how it is connected to a novel out-of-equilibrium statistical physics
measure that suppresses the confounding role of exponentially many
deep and isolated configurations (local minima of the error function)
and also amplifies the statistical weight of rare but extremely dense
regions of minima. We call this measure the *Robust Ensemble*
(RE). Moreover, we show that the RE allows us to derive novel and
exceptionally effective algorithms. One of these algorithms is closely
related to a recently proposed stochastic learning protocol used in
complex deep artificial neural networks NIPS2015_5761 , implying
that the underlying geometrical structure of the RE may provide an
explanation for its effectiveness.

In the present work, we consider discrete neural networks with only one or two layers, which can be studied analytically. However, we believe that these results should extend to deep neural networks of which the models studied here are building blocks, and in fact to other learning problems as well. We are currently investigating this work-in-progress .

## Ii Interacting replicas as a tool for seeking dense regions

In statistical physics, the canonical ensemble describes the equilibrium (i.e., long-time limit) properties of a stochastic process in terms of a probability distribution over the configurations

of the system , where is the energy of the configuration, an inverse temperature, and the normalization factor is called the partition function and can be used to derive all thermodynamic properties. This distribution is thus defined whenever a function is provided, and indeed it can be studied and provide insight even when the system under consideration is not a physical system. In particular, it can be used to describe interesting properties of optimization problems, in which has the role of a cost function that one wishes to minimize; in these cases, one is interested in the limit , which corresponds to assigning a uniform weight over the global minima of the energy function. This kind of description is at the core of the well-known Simulated Annealing algorithm kirkpatrick1983optimization .In the past few decades, equilibrium statistical physics descriptions have emerged as fundamental frameworks for studying the properties of a variety of systems which were previously squarely in the domain of other disciplines. For example, the study of the phase transitions of the random

-satisfiability problem (SAT) was linked to the algorithmic difficulty of finding solutions mezard2002analytic ; krzakala-csp . It was shown that the system can exhibit different phases, characterized by the arrangement of the space of solutions in one, many or a few connected clusters. Efficient (polynomial-time) algorithms appear to exist only if the system has so-called “unfrozen” clusters: extensive and connected regions of solutions in which at most a sub-linear number of variables have a fixed value. If, on the contrary, all solutions are “frozen” (belonging to clusters in which an extensive fraction of the variables take a fixed value, thus confined to sub-spaces of the space of configurations), no efficient algorithms are known. In the limiting case in which all variables—except at most a sub-linear number of them—are fixed, the clusters are isolated point-like regions, and the solutions are called “locked”.zdeborova2008lockedFor learning problems with discrete synapses, numerical experiments indicate that efficient algorithms also seek unfrozen solutions. In ref.

baldassi_subdominant_2015 , we showed that the equilibrium description in these problems is insufficient, in the sense that it predicts that the problem is always in the completely frozen phase in which all solutions are locked huang2014origin , in spite of the fact that efficient algorithms seem to exist. This motivated us to introduce a different measure, which ignores isolated solutions and enhances the statistical weight of large, accessible regions of solutions:(1) |

Here is a parameter that has the formal role of an inverse temperature and is a “local free entropy”:

(2) |

where some monotonically increasing function
of the distance between configurations, defined according to the model
under consideration. In the limit , this expression
reduces (up to an additive constant) to a “local entropy”: it
counts the number of minima of the energy, weighting them (via the
parameter ) by the distance from a reference configuration
. Therefore, if is large, only the configurations
that are surrounded by an exponential number of local minima will
have a non-negligible weight. By increasing the value of ,
it is possible to focus on narrower neighborhoods around ,
and at large values of the reference will also
with high probability share the properties of the surrounding minima.
This is illustrated in Fig. 1.
These large-deviation statistics seem to capture very well the behavior
of efficient algorithms on discrete neural networks, which invariably
find solutions belonging to high-density regions when these regions
exist, and fail otherwise. These solutions therefore could be *rare*
(i.e., not emerge in a standard equilibrium description), and yet
be *accessible* (i.e., there exist efficient algorithms that
are able to find them), and they are inherently *robust* (they
are immersed in regions of other “good” configurations). As discussed
in baldassi_subdominant_2015 , there is a relation between
the robustness of solutions in this sense and their good generalization
ability: this is intuitively understood in a Bayesian framework by
considering that a robust solution acts as a representative of a whole
extensive region.

It is therefore natural to consider using our large-deviation statistics as a starting point to design new algorithms, in much the same way that Simulated Annealing uses equilibrium statistics. Indeed, this was shown to work well in ref. baldassi_local_2016

. The main difficulty of that approach was the need to estimate the local (free) entropy

, which was addressed there using the Belief Propagation (BP) algorithm mezard_information_2009 .Here we demonstrate an alternative, general and much simpler approach. The key observation is that, when is a non-negative integer, we can rewrite the partition function of the large deviation distribution eq. (1) as:

This partition function describes a system of interacting replicas of the system, one of which acts as reference while the remaining are identical, subject to the energy and the interaction with the reference . Studying the equilibrium statistics of this system and tracing out the replicas is equivalent to studying the original large deviations model. This provides us with a very simple and general scheme to direct algorithms to explore robust, accessible regions of the energy landscape: replicating the model, adding an interaction term with a reference configuration and running the algorithm over the resulting extended system.

In fact, in most cases, we can further improve on this scheme by tracing
out the reference instead, which leaves us with a system of identical
interacting replicas describing what we call the *robust ensemble*
(RE):

(4) | |||||

(5) |

In the following, we will demonstrate how this simple procedure can be applied to a variety of different algorithms: Simulated Annealing, Stochastic Gradient Descent, and Belief Propagation. In order to demonstrate the utility of the method, we will focus on the problem of training neural network models.

## Iii Neural network models

Throughout this paper, we will consider for simplicity one main kind of neural network model, composed of identical threshold units arranged in a feed-forward architecture. Each unit has many input channels and one output, and is parameterized by a vector of “synaptic weights”

. The output of each unit is given by where is the vector of inputs.Since we are interested in connecting with analytical results, for the sake of simplicity all our tests have been performed using binary weights, , where denotes a hidden unit and an input channel. We should however mention that all the results generalize to the case of weights with multiple bits of precision baldassi2016learning . We denote by

the total number of synaptic weights in the network, which for simplicity is assumed to be odd. We studied the random classification problem: given a set of

random input patterns , each of which has a corresponding desired output , we want to find a set of parameters such that the network output equals for all patterns. Thus, for a single-layer network (also known as a perceptron), the condition could be written as

, where if and otherwise. For a fully-connected two-layer neural network (also known as committee or consensus machine), the condition could be written as (note that this assumes that all weights in the output unit are , since they are redundant in the case of binary ’s). A three-layer fully connected network would need to satisfy , and so on. In this work, we limited our tests to one- and two-layer networks.In all tests, we extracted all inputs and outputs in from unbiased, identical and independent distributions.

In order to use methods like Simulated Annealing and Gradient Descent, we need to associate an energy or cost to every configuration of the system

. One natural choice is just to count the number of errors (mis-classified patterns), but this is not a good choice for local search algorithms since it hides the information about what direction to move towards in case of error, except near the threshold. Better results can be obtained by using the following general definition instead: we define the energy

associated to each pattern as the minimum number of synapses that need to be switched in order to classify the pattern correctly. The total energy is then given by the sum of the energy for each pattern, . In the single layer case, the energy of a pattern is thus , where . Despite the simple definition, the expression for the two-layer case is more involved and is provided in the Appendix A.3. For deeper networks, the definition is conceptually too naive and computationally too hard, and it should be amended to reflect the fact that more than one layer is affected by training, but this is beyond the scope of the present work.We also need to define a distance function between replicas of the system. In all our tests, we used the squared distance , which is proportional to the Hamming distance in the binary case.

## Iv Replicated Simulated Annealing

We claim that there is a general strategy which can be used by a system of interacting replicas to seek dense regions of its configuration space. The simplest example of this is by sampling the configuration space with a Monte Carlo method moore2011nature which uses the objective functions given by eqs. (II) or (4), and lowering the temperature via a Simulated Annealing (SA) procedure, until either a zero of the energy (“cost”) or a “give-up condition” is reached. For simplicity, we use the RE, in which the reference configuration is traced out (eq. (4)), and we compare our method to the case in which the interaction between the replicas is absent (i.e. , which is equivalent to running parallel independent standard Simulated Annealing algorithms on the cost function). Besides the annealing procedure, in which the inverse temperature is gradually increased during the simulation, we also use a “scoping” procedure, which consists in gradually increasing the interaction , with the effect of reducing the average distance between the replicas. Intuitively, this corresponds to exploring the energy landscape on progressively finer scales (Fig. 1).

Additionally, we find that the effect of the interaction among replicas
can be almost entirely accounted for by adding a prior on the choice
of the moves within an otherwise standard Metropolis scheme, while
still maintaining the detailed balance condition (of course, this
reduces to the standard Metropolis rule for ). The elementary
Monte Carlo iteration step can therefore be summarized as follows:
*i)* we choose a replica uniformly at random; *ii)* we choose
a candidate weight to flip for that replica according to a prior computed
from the state of all the replicas; *iii)* we estimate the single-replica
energy cost of the flip and accept it according to a standard Metropolis
rule. The sampling technique (step *ii* above) and the parameters
used for the simulations are described in the Appendix B.

In Fig. 2, we show the results for the perceptron; an analogous figure for the committee machine, with similar results, is shown in the Appendix, Fig. 7. The analysis of the scaling with demonstrates that the interaction is crucial to finding a solution in polynomial time: the non-interacting version scales exponentially and it rapidly becomes impossible to find solutions in reasonable times. Our tests also indicate that the difference in performance between the interacting and non-interacting cases widens greatly with increasing the number of patterns per synapse . As mentioned above, this scheme bears strong similarities to the Entropy-driven Monte Carlo (EdMC) algorithm that we proposed in baldassi_local_2016 , which uses the Belief Propagation algorithm (BP) to estimate the local entropy around a given configuration. The main advantage of using a replicated system is that it avoids the need to use BP, which makes the procedure much simpler and more general. On the other hand, in systems where BP is able to provide a reasonable estimate of the local entropy, it can do so directly at a given temperature, and thus avoids the need to thermalize the replicas. Therefore, the landscapes explored by the replicated SA and EdMC are in principle different, and it is possible that the latter has fewer local minima; this however does not seem to be an issue for the neural network systems considered here.

## V Replicated Gradient Descent

Monte Carlo methods are computationally expensive, and may be infeasible for large systems. One simple alternative general method for finding minima of the energy is using Gradient Descent (GD) or one of its many variants. All these algorithms are generically called backpropagation algorithms in the neural networks context

rumelhart1988learning. Indeed, GD—in particular Stochastic GD (SGD)—is the basis of virtually all recent successful “deep learning” techniques in Machine Learning. The two main issues with using GD are that it does not in general offer any guarantee to find a global minimum, and that convergence may be slow (in particular for some of the variables, cf. the “vanishing gradient” problem

vanishing-gradient which affects deep neural network architectures). Additionally, when training a neural network for the purpose of inferring (generalizing) a rule from a set of examples, it is in general unclear how the properties of the local minima of the energy on the training set are related to the energy of the test set, i.e., to the generalization error.GD is defined on differentiable systems, and thus it cannot be applied directly to the case of systems with discrete variables considered here. One possible work-around is to introduce a two-level scheme, consisting in using two sets of variables, a continuous one and a discrete one , related by a discretization procedure , and in computing the gradient over the discrete set but adding it to the continuous set: (where is a gradient step, also called learning rate in the machine learning context). For the single-layer perceptron with binary synapses, using the energy definition provided above, in the case when the gradient is computed one pattern at a time (in machine learning parlance: using SGD with a minibatch size of ), this procedure leads to the so-called “Clipped Perceptron” algorithm (CP). This algorithm is not able to find a solution to the training problem in the case of random patterns, but simple (although non-trivial) variants of it are (SBPI and CP+R, see baldassi-et-all-pnas ; baldassi-2009 ). In particular CP+R was adapted to two-layer networks (using a a simplified version of the two-level SGD procedure described above) and was shown in baldassi_subdominant_2015

to be able to achieve near-state-of-the-art performance on the MNIST database

lecun1998gradient . The two-level SGD approach was also more recently applied to multi-layer binary networks with excellent results in courbariaux2015binaryconnect ; courbariaux2016binarynet , along with an array of additional heuristic modifications of the SGD algorithm that have become standard in application-driven works (e.g., batch renormalization). In those cases, however, the back-propagation of the gradient was performed differently, either because the output of each unit was not binary courbariaux2015binaryconnect or as a work-around for the use of a different definition for the energy, which required the introduction of additional heuristic mechanisms courbariaux2016binarynet .Almost all the above-mentioned results are purely heuristic (except in the on-line generalization setting, which is not considered in the present work). Indeed, even just using the two-level SGD is heuristic in this context. Nevertheless, here we demonstrate that, as in the case of Simulated Annealing of the previous section, replicating the system and adding a time-dependent interaction term, i.e., performing the gradient descent over the robust ensemble (RE) energy defined in eq. (5), leads to a noticeable improvement in the performance of the algorithm, and that when a solution is found it is indeed part of a dense region, as expected. We showed in baldassi_subdominant_2015 that solutions belonging to maximally dense regions have better generalization properties than other solutions; in other words, they are less prone to overfitting.

It is also important to note here that the stochastic nature of the SGD algorithm is essential in this case in providing an entropic component which counters the purely attractive interaction between replicas introduced by eq. (5), since the absolute minima of the replicated energy of eq. (4) are always obtained by collapsing all the replicas in the minima of the original energy function. Indeed, the existence of an entropic component allowed us to derive the RE definition by using the interaction strength to control the distance via a Legendre transform in the first place in eq. (2); in order to use this technique with a non-stochastic minimization algorithm, the distance should be controlled explicitly instead.

In Fig. 3 we show the results for a fully
connected committee machine,^{1}^{1}1The code used for these tests is available at codeRSGD .
demonstrating that the introduction of the interaction term greatly
improves the capacity of the network (from to almost
patterns per synapse), finds configurations with a lower error rate
even when it fails to solve the problem, and generally requires fewer
presentations of the dataset (epochs). The graphs show the results
for replicas in which the gradient is computed for every
patterns (the so-called minibatch size); we observed the same general
trend for all cases, even with minibatch sizes of (in the Appendix,
Fig. 8, we show the results for
and minibatch size ). We also observed the same effect in the
perceptron, although with less extensive tests, where this algorithm
has a capacity of at least . All technical details are provided
in Appendix C. These results are in perfect agreement
with the analysis of the next section, on Belief Propagation, which
suggests that this replicated SGD algorithm has near-optimal capacity.

It is interesting to note that a very similar approach—a replicated system in which each replica is attracted towards a reference configuration, called Elastic Averaged SGD (EASGD)—was proposed in NIPS2015_5761 (see also zhang2016distributed ) using deep convolutional networks with continuous variables, as a heuristic way to exploit parallel computing environments under communication constraints. Although it is difficult in that case to fully disentangle the effect of replicating the system from the other heuristics (in particular the use of “momentum” in the GD update), their results clearly demonstrate a benefit of introducing the replicas in terms of training error, test error and convergence time. It seems therefore plausible that, despite the great difference in complexity between their network and the simple models studied in this paper, the general underlying reason for the effectiveness of the method is the same, i.e., the existence of accessible robust low-energy states in the space of configurations work-in-progress .

## Vi Replicated Belief Propagation

Belief Propagation (BP), also known as Sum-Product, is an iterative message-passing method that can be used to describe a probability distribution over an instance described by a factor graph in the correlation decay approximation mackay2003information ; yedidia2005constructing . The accuracy of the approximation relies on the assumption that, when removing an interaction from the network, the nodes involved in that interaction become effectively independent, an assumption linked to so-called Replica Symmetry (RS) in statistical physics.

Algorithmically, the method can be briefly summarized as follows.
The goal is to solve a system of equations involving quantities called
messages. These messages represent single-variable *cavity* marginal
probabilities, where the term “cavity” refers to the fact that
these probabilities neglect the existence of a node in the graph.
For each edge of the graph there are two messages going in opposite
directions. Each equation in the system gives the expression of one
of the messages as a function of its neighboring messages. The expressions
for these functions are derived from the form of the interactions
between the variables of the graph. The resulting system of equations
is then solved iteratively, by initializing the messages in some arbitrary
configuration and updating them according to the equations until a
fixed point is eventually reached. If convergence is achieved, the
messages can be used to compute various quantities of interest; among
those, the most relevant in our case are the single-site marginals
and the thermodynamic quantities such as the entropy and the free
energy of the system. From single-site marginals, it is also possible
to extract a configuration by simply taking the
of each marginal.

In the case of our learning problem, the variables are the synaptic weights, and each pattern represents a constraint (i.e., an “interaction”) between all variables, and the form of these interactions is such that BP messages updates can be computed efficiently. Also note that, since our variables are binary, each of the messages and marginals can be described with a single quantity: we generally use the magnetization, i.e., the difference between the probability associated to the states and . We thus generally speak of the messages as “cavity magnetizations” and the marginals as “total magnetizations”. The factor graph, the BP equations and the procedure to perform the updates efficiently are described in detail in the Appendix section D.1, closely following the work of braunstein-zecchina . Here, we give a short summary. We use the notation for the message going from the node representing the weight variable to the node representing pattern at iteration step , and for the message going in the opposite direction, related by the BP equation:

(6) |

where indicates the set of all factor nodes connected to . The expressions for the total magnetizations are identical except that the summation index runs over the whole set . A configuration of the weights is obtained from the total magnetizations simply as . The expression for as a function of is more involved and it’s reported in the Appendix section D.1.

As mentioned above, BP is an inference algorithm: when it converges, it describes a probability distribution. One particularly effective scheme to turn BP into a solver is the addition of a “reinforcement” term braunstein-zecchina : a time-dependent local field is introduced for each variable, proportional to its own marginal probability as computed in the previous iteration step, and is gradually increased until the whole system is completely biased toward one configuration. This admittedly heuristic scheme is quite general, leads to very good results in a variety of different problems, and can even be used in cases in which unmodified BP would not converge or would provide a very poor approximation (see e.g. bailly2011finding ). In the case of the single layer binary network such as those considered in this paper, it can reach a capacity of patterns per synapse braunstein-zecchina , which is consistent with the value at which the structure of solution-dense regions breaks baldassi_subdominant_2015 .

The reason for the effectiveness of reinforced BP has not been clear. Intuitively, the process progressively focuses on smaller and smaller regions of the configuration space, with these regions determined from the current estimate of the distribution by looking in the “most promising” direction. This process has thus some qualitative similarities with the search for dense regions described in the previous sections. This analogy can be made precise by writing the BP equations for the system described by eq. (II). There are in this case two equivalent approaches: the first is to use the local entropy as the energy function, using a second-level BP to estimate the local entropy itself. This approach is very similar to the so called 1-step replica-symmetry-breaking (1RSB) cavity equations (see mezard_information_2009 for a general introduction). The second approach is to replicate the system, considering vector variables of length , and assuming an internal symmetry for each variable, i.e. that all marginals are invariant under permutation of the replica indices, similarly to what is done in kabashima2005 : . The result in the two cases is the same (this will be shown in a technical work, in preparation, where the connection between the large deviations measure and the 1RSB equilibrium description is also made explicit). Since BP assumes replica symmetry, the resulting message passing algorithm reproduces quite accurately the analytical results at the RS level. As explained in baldassi_subdominant_2015 , these results can however become wrong, in particular for high values of , and , due to the onset of correlations (the so called replica-symmetry-breaking – RSB – effect mezard_information_2009 ). More specifically, in this model the RS solution assumes that there is a single dense region comprising the robust ensemble (RE), while the occurrence of RSB implies that there are several maximally dense regions. As a consequence this algorithm is not a very good candidate as a solver. A more correct description—which could then lead to a more controlled solver—would thus require a third level of BP equations, or equivalently an assumption of symmetry-breaking in the structure of the marginals .

Fortunately, it turns out that that there is a different way of applying BP to the replicated system, leading to an efficient solver which is both very similar to the reinforced BP algorithm and reasonably well described by the theoretical results. Instead of considering the joint distribution over all replicated variables at a certain site

, we simply replicate the original factor graph times; then, for each site , we add an extra variable , and interactions, between each variable and . Finally, since the problem is symmetric, we assume that each replica of the system behaves in exactly the same way, and therefore that the same messages are exchanged along the edges of the graph regardless of the replica index. This assumption neglects some of the correlations between the replicated variables, but allows us to work only with a single system, which is identical to the original one except that each variable now also exchanges messages with identical copies of itself through an auxiliary variable (which we can just trace away at this point, as in eq. (4)). The procedure is shown graphically in Fig. 4. At each iteration step , each variable receives an extra message of the form:(7) |

where is the cavity magnetization resulting from the rest of the factor graph at time , and is the interaction strength. This message simply enters as an additional term in the sum in the right-hand side of eq. (6). Note that, even though we started from a system of replicas, after the transformation we are no longer constrained to keep in the integer domain. The reinforced BP braunstein-zecchina , in contrast, would have a term of the form:

(8) |

The latter equation uses a single parameter instead
of two, and is expressed in terms of the total magnetization
instead of the cavity magnetization . Despite
these differences, these two terms induce exactly the same BP fixed
points if we set and
(see the Appendix section D.2). This result
is somewhat incidental, since in the context of our analysis the limit
should be reached gradually, otherwise the results
should become trivial (see figure 1),
and the reason why this does not happen when setting
in eq. (7) is only related to the approximations
introduced by the simplified scheme of figure (4).
As it turns out, however, choosing slightly different mappings with
both and diverging gradually (e.g.
and with increasing from to
) can lead to update rules with the same qualitative behavior
and very similar quantitative effects, such that the performances
of the resulting algorithms are hardly distinguishable, and such that
the aforementioned approximations do not thwart the consistency with
the theoretical analysis. This is shown and discussed Appendix D.2.
Using a protocol in which is increased gradually, rather
than being set directly at , also allows the algorithm to
reach a fixed point of the BP iterative scheme before proceeding to
the following step, which offers more control in terms of the comparison
with analytical results, as discussed in the next paragraph. In this
sense, we therefore have derived a qualitative explanation of the
effectiveness of reinforced BP within a more general scheme for the
search of accessible dense states. We call this algorithm Focusing
BP (fBP).^{2}^{2}2The code implementing fBP on Committee Machines with binary synaptic
weights is available at codeFBP .

Apart from the possibility of using fBP as a solver, by gradually increasing and until a solution is found, it is also interesting to compare its results at fixed values of and with the analytical predictions for the perceptron case which were derived in baldassi_subdominant_2015 ; baldassi_local_2016 , since at the fixed point we can use the fBP messages to compute such values as the local entropy, the distance from the reference configuration, and the average overlap between replicas (defined as for any and , where denotes the thermal averaging). The expressions for these quantities are provided in Appendix D.1. The results indicate that fBP is better than the alternatives described above in overcoming the issues arising from RSB effects. The 1RSB scheme describes a non-ergodic situation which arises from the break-up of the space of configurations into a collection of several separate equivalent ergodic clusters, each representing a possible state of the system. These states are characterized by the typical overlap inside the same state (the intra-cluster overlap ) and the typical overlap between configurations from any two different states (the inter-cluster overlap ). Fig. 5 shows that the average overlap between replicas computed by the fBP algorithm transitions from to when is increased (and the distance with the reference is decreased as a result). This suggests that the algorithm has spontaneously chosen one of the possible states of high local entropy in the RE, achieving an effect akin to the spontaneous symmetry breaking of the 1RSB description. Within the state, replica symmetry holds, so that the algorithm is able to eventually find a solution to the problem. Furthermore, the resulting estimate of the local entropy is in very good agreement with the 1RSB predictions up to at least (see the Appendix, figure 11).

Therefore, although this algorithm is not fully understood from the theoretical point of view, it offers a valuable insight into the reason for the effectiveness of adding a reinforcement term to the BP equations. It is interesting in this context to observe that the existence of a link between the reinforcement term and the robust ensemble seems consistent with some recent results on random bicoloring constraint satisfaction problems braunstein_large_2016 , which showed that reinforced BP finds solutions with shorter “whitening times” with respect to other solvers: this could be interpreted as meaning they belong to a region of higher solution density, or are more central in the cluster they belong to. Furthermore, our algorithm can be used to estimate the point up to which accessible dense states exist, even in cases, like multi-layer networks, where analytical calculations are prohibitively complex.

Fig. 6 shows the result of experiments performed on a committee machine with the same architecture and same of Fig. 3. The implementation closely follows braunstein-zecchina with the addition of the self-interaction eq. (7), except that great care is required to correctly estimate the local entropy at large , due to numerical issues (see Appendix D.1). The figure shows that fBP finds that dense states (where the local entropy curves approach the upper bound at small distances) exist up to nearly patterns per synapse, and that when it finds those dense states it is correspondingly able to find a solution, in perfect agreement with the results of the replicated gradient descent algorithm.

Finally, we also performed some exploratory tests applying fBP on the random -satisfiability problem, and we found clear indications that the performance of this algorithm is similar to that of Survey-Propagation-based algorithms mezard2002analytic ; marino2015backtrackingSP , although a more detailed analysis is required to draw a more precise comparison.. The results are reported in Appendix D.4.

## Vii Discussion

In this paper, we have presented a general scheme that can be used to bias the search for low-energy configurations, enhancing the statistical weight of large, accessible states. Although the underlying theoretical description is based on a non-trivial large deviation measure, its concrete implementation is very simple—replicate the system and introduce an interaction between the replicas—and versatile, in that it can be generally applied to a number of different optimization algorithms or stochastic processes. We demonstrated this by applying the method to Simulated Annealing, Gradient Descent and Belief Propagation, but it is clear that the list of possible applications may be much longer. The intuitive interpretation of the method is also quite straightforward: a set of coupled systems is less likely to get trapped in narrow minima, and will instead be attracted to wide regions of good (and mostly equivalent) configurations, thus naturally implementing a kind of robustness to details of the configurations.

The utility of this kind of search depends on the details of the problem under study. Here we have mainly focused on the problem of training neural networks, for a number of reasons. The first is that, at least in the case of single-layer networks with discrete synapses, we had analytical and numerical evidence that dense, accessible states exist and are crucial for learning and improving the generalization performance, and we could compare our findings with analytical results. The second is that the general problem of training neural networks has been addressed in recent years via a sort of collective search in the space of heuristics, fueled by impressive results in practical applications and mainly guided by intuition; heuristics are evaluated based on their effectiveness in finding accessible states with good generalization properties. It seems reasonable to describe these accessible states as regions of high local entropy, i.e., wide, very robust energy minima: the center of such a region can act as a Bayesian estimator for the whole extensive neighborhood. Here we showed a simple way to exploit the existence of such states efficiently, whatever the optimization algorithm used. This not only sheds light on previously known algorithms, but also suggests improvements or even entirely new algorithms. Further work is required to determine whether the same type of phenomenon that we observed here in simple models actually generalizes to the deep and complex networks currently used in machine learning applications (the performance boost obtained by the EASGD algorithm of NIPS2015_5761 being a first indication in this direction), and to investigate further ways to improve the performance of learning algorithms, or to overcome constraints (such as being limited to very low-precision computations).

It is also natural to consider other classes of problems in which this analysis may be relevant. One application would be solving other constraint satisfaction problems. For example, in baldassi_local_2016 we demonstrated that the EdMC algorithm can be successfully applied to the random -satisfiability problem, even though we had to resort to a rough estimate of the local entropy due to replica symmetry breaking effects. We have shown here that the fBP algorithm presented above is also effective and efficient on the same problem. It is also interesting to note here that large-deviation analyses—although different from the one of the present paper—of a similar problem, the random bicoloring constraint satisfaction problem,dall2008entropy ; braunstein_large_2016 have shown that atypical “unfrozen” solutions exist (and can be found with the reinforced BP algorithm) well beyond the point in the phase diagram where the overwhelming majority of solutions have become “frozen”. Finally, an intriguing problem is the development of a general scheme for a class of out-of-equilibrium processes attracted to accessible states: even when describing a system which is unable to reach equilibrium in the usual thermodynamic sense or is driven by some stochastic perturbation, it is still likely that its stationary state can be characterized by a large local entropy.

###### Acknowledgements.

We wish to thank Y. LeCun and L. Bottou for encouragement and interesting discussions about future directions for this work. CBa, CL and RZ acknowledge the European Research Council for grant n° 267915.## Appendix A Model and notation

This Appendix text contains all the technical details of the algorithms described in the main text, the techniques and the parameters we used to obtain the results we reported. We also report some additional results and report other minor technical considerations.

Preliminarily, we set a notation used throughout the rest of this document which is slightly different from the one of the main text, but more suitable for this technical description.

### a.1 The network model

As described in the main text, we consider an ensemble of neural networks with

units and binary variables

where is the unit index, is the synaptic index and is the replica index. Each network has thus synapses, where is divisible by . For simplicity, we assume both and to be odd. The output of each unit is defined by a function . The output of the network is defined by a function where represents the input to the -th unit. In the case , this is equivalent to a single-layer network (also known as perceptron). In the case where all are identical for each , this is equivalent to a fully-connected two-layer network (also known as committee machine or consensus machine). If the are different for different values of , this is a tree-like committee machine. Note that, due to the binary constraint on the model, adding weights to the second layer is redundant, since for all negative weights in the second layer we could always flip both its weight and all the weights of the unit connected to it. Therefore, without loss of generality, we just set the weights of the second layer to , resulting in the above definition of the output function .The scalar product between two replicas and is defined as . For brevity of notation, in cases where the unit index does not play a role, we will often just use a single index , e.g. .

### a.2 Patterns

The networks are trained on random input/output associations, i.e. patterns, where is the pattern index. The parameter determines the load of the network, so that the number of patterns is proportional to the number of synapses. The inputs are binary vectors of elements with entries , and the desired outputs are also binary, . Both the inputs and the outputs are extracted at random and are independent and identically distributed (i.i.d.), except in the case of the fully-connected committee machine where for all and therefore we only extract the values for .

We also actually exploit a symmetry in the problem and set all desired outputs to , since for each pattern its opposite must have an opposite output, i.e. we can always transform an input output pair into , where the new pattern has the same probability as .

### a.3 Energy definition

The energy, or cost, for each pattern is defined as the minimum number of weights which need to be switched in order to correctly classify the pattern, i.e. in order to satisfy the relation . The total energy is the sum of the energies for all patterns, .

If the current configuration of the weights satisfies the pattern, the corresponding energy is obviously . Therefore, if the training problem is satisfiable, the ground states with this energy definition are the same as for the easier energy function given in terms of the number of errors.

If the current configuration violates the pattern, the energy can be computed as follows: we need to compute the minimum number of units of the first level which need to change their outputs, choose the units which are easiest to fix, and for each of them compute the minimum number of weights which need to be changed. In formulas:

(9) |

where:

(10) | |||||

(11) | |||||

(12) | |||||

(13) |

where the function returns its argument sorted in ascending order. The above auxiliary quantities all depend on , but we omitted the dependency for clarity of notation.

In the single-layer case the expression simplifies considerably, since and reduces to .

## Appendix B Replicated Simulated Annealing

We run Simulated Annealing (plus “scoping”) on a system of interacting replicas. For simplicity, we trace away the reference configuration which mediates the interaction. Thus, at any given step, we want to sample from a probability distribution

(14) | |||||

The reference configuration is traced out in this representation, but we can obtain its most probable value by just computing . It is often the case that, when the parameters are chosen appropriately, , i.e. that the energy of the center is lower than that of the group of replicas. In fact, we found this to be a good rule-of-thumb criterion to evaluate the choice of the parameters in the early stages of the algorithmic process.

The most straightforward way to perform the sampling (at fixed and ) is by using the Metropolis rule; the proposed move is to flip a random synaptic weight from a random replica. Of course the variation of the energy associated to the candidate move now includes the interaction term, parametrized by , which introduces a bias that favors movements in the direction of the center of mass of the replicas.

We also developed an alternative rule for choosing the moves in a biased way which implicitly accounts for the interaction term while still obeying the detailed balance condition. This alternative rule is generally valid in the presence of an external field and is detailed at the end of this section. Its advantage consists in reducing the rejection rate, but since the move proposal itself becomes more time consuming it is best suited to systems in which computing the energy cost of a move is expensive, so its usefulness depends on the details of the model.

### b.1 Computing the energy shifts efficiently

Here we show how to compute efficiently the quantity when and only differ in the value of one synaptic weight and the energy is defined as in eq. (9). To this end, we define some auxiliary quantities in addition to the ones required for the energy computation, eqs. (10)-(13) (note that we omit the replica index here since this needs to be done for each replica independently):

(15) | |||||

(16) | |||||

(17) |

These quantities must be recomputed each time a move is accepted, along with (10)-(13). Note however that in later stages of the annealing process most moves are rejected, and the energy shifts can be computed very efficiently as we shall see below.

Preliminarily, we note that any single-flip move only affects the energy contribution from patterns in .

The contribution to the energy shift for a proposed move is most easily written in pseudo-code:

Indeed, this function is greatly simplified in the single-layer case .

### b.2 Efficient Monte Carlo sampling

Here we describe a Monte Carlo sampling method which is a modification of the Metropolis rule when the system uses binary variables and the Hamiltonian function can be written as:

(18) |

where the external fields can only assume a finite (and much smaller than ) set of values. The factor is introduced merely for convenience. Comparing this to eq. (14), we see that, having chosen a replica index uniformly at random, we can identify

(19) |

.

Given a transition probability to go from state to state , , the detailed balance equation reads:

(20) |

Let us split the transition explicitly in two steps: choosing the index and accepting the move. The standard Metropolis rule is: pick an index uniformly at random, propose the flip of , accept it with probability , where . We want to reduce the rejection rate and incorporate the effect of the field in the proposal instead. We write:

(21) |

where is the choice of the index, and is the acceptance of the move. Usually is uniform and we ignore it, but here instead we try to use it to absorb the external field term in the probability distribution. From detailed balance we have:

(22) | |||||

so if we could satisfy:

(23) |

then the acceptance would simplify to the usual Metropolis rule, involving only the energy shift . This will turn out to be impossible, yet easily fixable, so we still first derive the condition implied by eq. (23). The key observation is that there is a finite number of classes of indices in , based on the limited number of values that can take (in the case of eq. (19) there are possible values). Let us call the possible classes, such that , and let us call their sizes, with the normalization condition that . Within a class, we must choose the move uniformly.

Then is determined by the probability of picking a class, which in principle could be a function of all the values of the : . Suppose now that we have picked an index in a class . The transition to would bring it into class , and the new class sizes would be

therefore:

(24) |

Since the only values of directly involved in this expression are and , it seems reasonable to restrict the dependence of and only on those values. Let us also call , which is unaffected by the transition. Therefore we can just write:

(25) |

Furthermore, we can assume – purely for simplicity – that:

(26) |

which allows us to restrict ourselves in the following to the case , and which implies that the choice of the index will proceed like this: we divide the indices in super-classes of size and we choose one of those according to their size; then we choose either the class or according to ; finally, we choose an index inside the class uniformly at random. Considering this process, what we actually need to determine is the conditional probability of choosing once we know we have chosen the super-class :

(27) |

Looking at eq. (23) we are thus left with the condition:

(28) |

Considering that we must have , this expression allows us to compute recursively for all values of . The computation can be carried out analytically and leads to where the function is defined as:

(29) |

with the hypergeometric function. However, we should also have , while and therefore this condition is only satisfied for (in which case we recover

, i.e. the standard uniform distribution, as expected).

Therefore, as anticipated, eq. (23) can not be satisfied^{3}^{3}3Strictly speaking we have not proven this, having made some assumptions
for simplicity. However it is easy to prove it in the special case
in which , since then our assumptions
become necessary., and we are left with a residual rejection rate for the case .
This is reasonable, since in the limit of very large (i.e. very
large in the case of eq. (19)) the probability
distribution of each spin must be extremely peaked on the state in
which all replicas are aligned, such that the combined probability
of all other states is lower than the probability of staying in the
same configuration. Therefore we have (still for ):

(30) | |||||

(31) |

where is the Kronecker delta symbol. The last condition can be satisfied by choosing a general acceptance rule of this form:

(32) |

where

In practice, the effect of this correction is that the state where all the variables in class are already aligned in their preferred direction is a little “clingier” than the others, and introduces an additional rejection rate (which however is tiny when either is small or is large).

The final procedure is thus the following: we choose a super-class at random with probability , then we choose either or according to and finally pick another index uniformly at random within the class.

This procedure is highly effective at reducing the rejection rate induced by the external fields. As mentioned above, depending on the problem, if the computation of the energy shifts is particularly fast, it may still be convenient in terms of CPU time to produce values uniformly and rejecting many of them, rather then go through a more involved sampling procedure. Note however that the bookkeeping operations required for keeping track of the classes compositions and their updates can be performed efficiently, in time with space, by using an unsorted partition of the spin indices (which allows for efficient insertion/removal) and an associated lookup table. Therefore, the additional cost of this procedure is a constant factor at each iteration.

Also, the function involves the evaluation of a hypergeometric function, which can be relatively costly; its values however can be pre-computed and tabulated if the memory resources allow it, since they are independent from the problem instance. For large values of , it can also be efficiently approximated by a series expansion. It is convenient for that purpose to change variables to

(note that ). We give here for reference the expansion up to , which ensures a maximum error of for :

(33) |

Finally, note that the assumption of eq. (26)
is only justified by simplicity; it is likely that a different choice
could lead to a further improved dynamics.^{4}^{4}4We did in fact generalize and improve this scheme after the preparation
of this manuscript, see baldassi2016rrrmc .

### b.3 Numerical simulations details

Our Simulated Annealing procedure was performed as follows: we initialized the replicated system in a random configuration, with all replicas being initialized equally. The initial inverse temperature was set to , and the initial interaction strength to . We then ran the Monte Carlo simulation, choosing a replica index at random at each step and a synaptic index according to the modified Metropolis rule described in the previous section, increasing both and , by a factor and respectively, for each accepted moves. The gradual increase of is called ‘annealing’ while the gradual increase of is called ‘scoping’. Of course, since with our procedure the annealing/scoping step is fixed, the quantities and should scale with . The simulations are stopped as soon as any one of the replicas reaches zero energy, or after consecutive non-improving moves, where a move is classified as non-improving if it is rejected by the Metropolis rule or it does not lower the energy (this definition accounts for the situation where the system is trapped in a local minimum with neighboring equivalent configurations at large , in which case the algorithm would keep accepting moves with without doing anything useful).

In order to compare our method with standard Simulated Annealing, we just removed the interaction between replicas from the above described case, i.e. we set . This is therefore equivalent to running independent (except for the starting configurations) procedures in parallel, and stopping as soon as one of them reaches a solution.

In order to determine the scaling of the solution time with , we followed the following procedure: for each sample (i.e. patterns assignment) we ran the algorithm with different parameters and recorded the minimum number of iterations required to reach a solution. We systematically explored these values of the parameters: , , , (the latter two only in the interacting case, of course). This procedure gives us an estimate for the minimum number of iterations required to solve a typical problem at a given value of , and . We tested samples for each value of . Since the interacting case has additional parameters, this implies that there were more optimization opportunities, attributable to random chance; this however is not remotely sufficient to explain the difference in performance between the two cases: in fact, comparing instead for the typical value of iterations required (i.e. optimizing the average iterations over ) gives qualitatively similar results, since once a range of good values for the parameters is found the iterations required to reach a solution are rather stable across samples.

The results are shown in figure 2 of the main text for the single-layer case at and figure 7 for the fully-connected two-layer case (committee machine) at and . In both cases we used , which seems to provide good results (we did not systematically explore different values of ). The values of were chosen so that the standard SA procedure would be able to solve some instances at low in reasonable times (since the difference in performance between the interacting and non-interacting cases widens greatly with increasing ). The results show a different qualitative behavior in both cases, polynomial for the interacting case and exponential for the non-interacting cases. All fits were performed directly in logarithmic scale. A similar behavior is observed for the tree-like committee machine (not shown).

## Appendix C Replicated Gradient Descent

### c.1 Gradient computation

As mentioned in the main text, we perform a stochastic gradient descent on binary networks using the energy function of eq. (9) by using two sets of variables: a set of continuous variables

and the corresponding binarized variables

, related by . We use the binarized variables to compute the energy and the gradient, and apply the gradient to the continuous variables. In formulas, the quantities at time are related to those at time by:(34) | |||||

(35) |

where is a learning rate and is a set of pattern indices (a so-called minibatch). A particularly simple scenario can be obtained by considering a single layer network without replication (, ) and a fixed learning rate, and by computing the gradient one pattern at a time (). In that case, where and the gradient is . Since the relation (35) is scale-invariant, we can just set and obtain

(36) |

where now the auxiliary quantities are discretized: if they are initialized as odd integers, they remain odd integers throughout the learning process. This is the so-called “Clipped Perceptron” (CP) rule, which is the same as the Perceptron rule (“in case of error, update the weights in the direction of the pattern, otherwise do nothing”) except that the weights are clipped upon usage to make them binary. Notably, the CP rule by itself does not scale well with ; it is however possible to make it efficient (see baldassi-et-all-pnas ; baldassi-2009 ).

In the two-layer case () the computation of the gradient is more complicated; it is however simpler than the computation of the energy shift which was necessary for Simulated Annealing (Algorithm 1), since we only consider infinitesimal variations when computing the gradient. The resulting expression is:

(37) |

i.e. the gradient is non-zero only in case of error, and only for those units which contribute to the energy computation (which turn up in the first terms of the sorted vector , see eqs. (10)-(13)). Again, since this gradient can take only possible values, we could set and use discretized odd variables for the .

It is interesting to point out that a slight variation of this update rule in which only the first, least-wrong unit is affected, i.e. in which the condition is changed to , was used in baldassi_subdominant_2015 , giving good results on a real-world learning task when a slight modification analogous to the one of baldassi-2009 was added. Note that, in the later stages of learning, when the overall energy is low, it is very likely that , implying that the simplification used in baldassi_subdominant_2015 likely has a negligible effect. The simplified version, when used in the continuous case, also goes under the name of “least action” algorithm mitchison1989bounds .

Having computed the gradient of for each system, the extension to the replicated system is rather straightforward, since the energy (with the traced-out center) becomes (cf. eqs. (4) and (5) in the main text):

(38) |

and therefore the gradient just has an additional term:

(39) |

Note that the trace operation brings the parameter into account. Using as control parameter, the update equation (34) for a replica becomes (we omit the unit index for simplicity):