Year: 2021

Abstract:
Message-passing algorithms based on the Belief Propagation (BP) equations constitute a well-known distributed computational scheme. They yield exact marginals on tree-like graphical models and have also proven to be effective in many problems defined on loopy graphs, from inference to optimization, from signal processing to clustering. The BP-based schemes are fundamentally different from stochastic gradient descent (SGD), on which the current success of deep networks is based. In this paper, we present and adapt to mini-batch training on GPUs a family of BP-based message-passing algorithms with a reinforcement term that biases distributions towards locally entropic solutions. These algorithms are capable of training multi-layer neural networks with performance comparable to SGD heuristics in a diverse set of experiments on natural datasets including multi-class image classification and continual learning, while being capable of yielding improved performances on sparse networks. Furthermore, they allow to make approximate Bayesian predictions that have higher accuracy than point-wise ones.

heading: 1 INTRODUCTION
text: Belief Propagation is a method for computing marginals and entropies in probabilistic inference problems (Bethe, 1935; Peierls, 1936; Gallager, 1962; Pearl, 1982). These include optimization problems as well once they are written as zero temperature limit of a Gibbs distribution that uses the cost function as energy. Learning is one particular case, in which one wants to minimize a cost which is a data dependent loss function. These problems are generally intractable and message-passing techniques have been particularly successful at providing principled approximations through efficient distributed computations.
A particularly compact representation of inference/optimization problems that is used to build massage-passing algorithms is provided by factor graphs. A factor graph is a bipartite graph composed of variables nodes and factor nodes expressing the interactions among variables. Belief Propagation is exact for tree-like factor graphs (Yedidia et al., 2003)), where the Gibbs distribution is naturally factorized, whereas it is approximate for graphs with loops. Still, loopy BP is routinely used with success in many real world applications ranging from error correcting codes, vision, clustering, just to mention a few. In all these problems, loops are indeed present in the factor graph and yet the variables are weakly correlated at long range and BP gives good results. A field in which BP has a long history is the statistical physics of disordered systems where it is known as Cavity Method (Mézard et al., 1987). It has been used to study the typical properties of spin glass models which represent binary variables interacting through random interactions over a given graph. It is very well known that in spin glass models defined on complete graphs and in locally tree-like random graphs, which are both loopy, the weak correlation conditions between variables may hold and BP give asymptotic exact results (Mézard & Montanari, 2009). Here we will mostly focus on neural networks ±1 binary weights and sign activation functions, for which the messages and the marginals can be described simply by the difference between the probabilities associated with the +1 and -1 states, the so called magnetizations. The effectiveness of BP for deep learning has never been numerically tested in a systematic way, however there is clear evidence that the weak correlation decay condition does not hold and thus BP convergence and approximation quality is unpredictable.
In this paper we explore the effectiveness of a variant of BP that has shown excellent convergence properties in hard optimization problems and in non-convex shallow networks. It goes under the
name of focusing BP (fBP) and is based on a probability distribution, a likelihood, that focuses on highly entropic wide minima, neglecting the contribution to marginals from narrow minima even when they are the majority (and hence dominate the Gibbs distribution). This version of BP is thus expected to give good results only in models that have such wide entropic minima as part of their energy landscape. As discussed in (Baldassi et al., 2016a), a simple way to define fBP is to add a "reinforcement" term to the BP equations: an iteration-dependent local field is introduced for each variable, with an intensity proportional to its marginal probability computed in the previous iteration step. This field is gradually increased until the entire system becomes fully biased on a configuration. The first version of reinforced BP was introduced in (Braunstein & Zecchina, 2006) as a heuristic algorithm to solve the learning problem in shallow binary networks. Baldassi et al. (2016a) showed that this version of BP is a limiting case of fBP, i.e., BP equations written for a likelihood that uses the local entropy function instead of the error (energy) loss function. As discussed in depth in that study, one way to introduce a likelihood that focuses on highly entropic regions is to create y coupled replicas of the original system. fBP equations are obtained as BP equations for the replicated system. It turns out that the fBP equations are identical to the BP equations for the original system with the only addition of a self-reinforcing term in the message passing scheme. The fBP algorithm can be used as a solver by gradually increasing the effect of the reinforcement: one can control the size of the regions over which the fBP equations estimate the marginals by tuning the parameters that appear in the expression of the reinforcement, until the high entropy regions reduce to a single configuration. Interestingly, by keeping the size of the high entropy region fixed, the fBP fixed point allows one to estimate the marginals and entropy relative to the region.
In this work, we present and adapt to GPU computation a family of fBP inspired message passing algorithms that are capable of training multi-layer neural networks with generalization performance and computational speed comparable to SGD. This is the first work that shows that learning by message passing in deep neural networks 1) is possible and 2) is a viable alternative to SGD. Our version of fBP adds the reinforcement term at each mini-batch step in what we call the Posterioras-Prior (PasP) rule. Furthermore, using the message-passing algorithm not as a solver but as an estimator of marginals allows us to make locally Bayesian predictions, averaging the predictions over the approximate posterior. The resulting generalization error is significantly better than those of the solver, showing that, although approximate, the marginals of the weights estimated by messagepassing retain useful information. Consistently with the assumptions underlying fBP, we find that the solutions provided by the message passing algorithms belong to flat entropic regions of the loss landscape and have good performance in continual learning tasks and on sparse networks as well.
We also remark that our PasP update scheme is of independent interest and can be combined with different posterior approximation techniques.
The paper is structured as follows: in Sec. 2 we give a brief review of some related works. In Sec. 3 we provide a detailed description of the message-passing equations and of the high level structure of the algorithms. In Sec. 4 we compare the performance of the message passing algorithms versus SGD based approaches in different learning settings.

heading: 2 RELATED WORKS
text: The literature on message passing algorithms is extensive, we refer to Mézard & Montanari (2009) and Zdeborová & Krzakala (2016) for a general overview. More related to our work, multilayer message-passing algorithms have been developed in inference contexts (Manoel et al., 2017; Fletcher et al., 2018), where they have been shown to produce exact marginals under certain statistical assumptions on (unlearned) weight matrices.
The properties of message-passing for learning shallow neural networks have been extensively studied (see Baldassi et al. (2020) and reference therein). Barbier et al. (2019) rigorously show that message passing algorithms in generalized linear models perform asymptotically exact inference under some statistical assumptions. Dictionary learning and matrix factorization are harder problems closely related to deep network learning problems, in particular to the modelling of a single intermediate layer. They have been approached using message passing in Kabashima et al. (2016) and Parker et al. (2014), although the resulting predictions are found to be asymptotically inexact (Maillard et al., 2021). The same problem is faced by the message passing algorithm recently proposed for a multi-layer matrix factorization scenario (Zou et al., 2021). Unfortunately, our framework as well
doesn’t yield asymptotic exact predictions. Nonetheless, it gives a message passing heuristic that for the first time is able to train deep neural networks on natural datasets, therefore sets a reference for the algorithmic applications of this research line.
A few papers advocate the success of SGD to the geometrical structure (smoothness and flatness) of the loss landscape in neural networks (Baldassi et al., 2015; Chaudhari et al., 2017; Garipov et al., 2018; Li et al., 2018; Pittorino et al., 2021; Feng & Tu, 2021). These considerations do not depend on the particular form of the SGD dynamics and should extend also to other types of algorithms, although SGD is by far the most popular choice among NNs practitioners due to its simplicity, flexibility, speed, and generalization performance.
While our work focuses on message passing schemes, some of the ideas presented here, such as the PasP rule, can be combined with algorithms for Bayesian neural networks’ training (HernándezLobato & Adams, 2015; Wu et al., 2018). Recent work extends BP by combining it with graph neural networks (Kuck et al., 2020; Satorras & Welling, 2021). Finally, some work in computational neuroscience shows similarities to our approach (Rao, 2007).

heading: 3 LEARNING BY MESSAGE PASSING

heading: 3.1 POSTERIOR-AS-PRIOR UPDATES
text: We consider a multi-layer perceptron with L hidden neuron layers, having weight and bias parameters W = {W `, b`}L`=0. We allow for stochastic activations P `(x`+1|z`), where z` is the neuron’s preactivation vector for layer `, and P ` is assumed to be factorized over the neurons. If no stochasticity is present, P ` just encodes an element-wise activation function. The probability of output y given an input x is then given by
P (y |x,W) = ∫ dx1:L L∏ `=0 P `+1(x`+1 | W `x` + b`), (1)
where for convenience we defined x0 = x and xL+1 = y. In a Bayesian framework, given a training set D = {(xn, yn)}n and a prior distribution over the weights qθ(W) in some parametric family, the posterior distribution is given by
P (W |D, θ) ∝ ∏ n P (yn |xn,W) qθ(W). (2)
Here ∝ denotes equality up to a normalization factor. Using the posterior one can compute the Bayesian prediction for a new data-point x through the average P (y |x, D, θ) =∫ dW P (y |x,W)P (W |D, θ). Unfortunately, the posterior is generically intractable due to the hard-to-compute normalization factor. On the other hand, we are mainly interested in training a distribution that covers wide minima of the loss landscape that generalize well (Baldassi et al., 2016a) and in recovering pointwise estimators within these regions. The Bayesian modeling becomes an auxiliary tool to set the stage for the message passing algorithms seeking flat minima. We also need a formalism that allows for mini-batch training to speed-up the computation and deal with large datasets. Therefore, we devise an update scheme that we call Posterior-as-Prior (PasP), where we evolve the parameters θt of a distribution qθt(W) computed as an approximate mini-batch posterior, in such a way that the outcome of the previous iteration becomes the prior in the following step. In the PasP scheme, θt retains the memory of past observations. We also add an exponential factor ρ, that we typically set close to 1, tuning the forgetting rate and playing a role similar to the learning rate in SGD. Given a mini-batch (Xt,yt) sampled from the training set at time t and a scalar ρ > 0, the PasP update reads
qθt+1(W) ≈ [ P (W |yt,Xt, θt) ]ρ , (3)
where ≈ denotes approximate equality and we do not account for the normalization factor. A first approximation may be needed in the computation of the posterior, a second to project the approximate posterior onto the distribution manifold spanned by θ (Minka, 2001). In practice, we will consider factorized approximate posterior in an exponential family and priors qθ in the same family, although Eq. 3 generically allow for more refined approximations.
Notice that setting ρ = 1, the batch-size to 1, and taking a single pass over the dataset, we recover the Assumed Density Filtering algorithm (Minka, 2001). For large enough ρ (including ρ = 1), the iterations of qθt will concentrate on a pointwise estimator. This mechanism mimics the reinforcement heuristic commonly used to turn Belief Propagation into a solver for constrained satisfaction problems (Braunstein & Zecchina, 2006) and related to flat-minima discovery (see focusing-BP in Baldassi et al. (2016a)). A different prior updating mechanism which can be understood as empirical Bayes has been used in Baldassi et al. (2016b).

heading: 3.2 INNER MESSAGE PASSING LOOP
text: While the PasP rule takes care of the reinforcement heuristic across mini-batches, we compute the mini-batch posterior in Eq. 3 using message passing approaches derived from Belief Propagation. BP is an iterative scheme for computing marginals and entropies of statistical models Mézard & Montanari (2009). It is most conveniently expressed on factor graphs, that is bipartite graphs where the two sets of nodes are called variable nodes and factor nodes. They respectively represent the variables involved in the statistical model and their interactions. Message from factor nodes to variable nodes and viceversa are exchanged along the edges of the factor graph for a certain number of BP iterations or until a fixed point is reached.
The factor graph for P (W |Xt,yt, θt) can be derived from Eq. 2, with the following additional specifications. For simplicity, we will ignore the bias term in each layer. We assume factorized qθt(W), each factor parameterized by its first two moments. In what follows, we drop the PasP iteration index t. For each example (xn, yn) in the mini-batch, we introduce the auxiliary variables x`n, ` = 1, . . . , L, representing the layers’ activations. For each example, each neuron in the network contributes a factor node to the factor graph. The scalar components of the weight matrices and the activation vectors become variable nodes. This construction is presented in Appendix A, where we also derive the message update rules on the factor graph. The factor graph thus defined is extremely loopy and straightforward iteration of BP has convergence issues. Moreover, in presence of a homogeneous prior over the weights, the neuron permutation symmetry in each hidden layer induces a strongly attractive symmetric fixed point that hinders learning. We work around these issues by breaking the symmetry at time t = 0 with an inhomogeneous prior. In our experiments a little initial heterogeneity is sufficient to obtain specialized neurons at each following time step. Additionally, we do not require message passing convergence in the inner loop (see Algorithm 1) but perform one or a few iterations for each θ update. We also include an inertia term commonly called damping factor in the message updates (see B.2). As we shall discuss, these simple rules suffice to train deep networks by message passing.
For the inner loop we adapt to deep neural networks four different message passing algorithms, all of which are well known to the literature although derived in simpler settings: Belief Propagation (BP), BP-Inspired (BPI) message passing, mean-field (MF), and approximate message passing (AMP). The last three algorithms can be considered approximations of the first one. In the following paragraphs we will discuss their common traits, present the BP updates as an example, and refer to Appendix A for an in-depth exposition. For all algorithms, message updates can be divided in a forward pass and backward pass, as also done in (Fletcher et al., 2018) in a multi-layer inference setting. The BP algorithm is compactly reported in Algorithm 1.
Meaning of messages. All the messages involved in the message passing can be understood in terms of cavity marginals or full marginals (as mentioned in the introduction BP is also known as Cavity Method, see (Mézard & Montanari, 2009)). Of particular relevance are m`ki and σ ` ki, denoting the mean and variance of the weights W `ki. The quantities x̂ ` in and ∆ ` in instead denote the mean and variance of the i-th neuron’s activation in layer ` for a given input xn.
Scalar free energies. All message passing schemes are conveniently expressed in terms of two functions that correspond to the effective free energy (Zdeborová & Krzakala, 2016) of a single
neuron and of a single weight respectively :
ϕ`(B,A, ω, V ) = log ∫ dxdz e− 1 2Ax 2+Bx P ` (x|z) e− (ω−z)2 2V ` = 1, . . . , L (4)
ψ(H,G, θ) = log ∫ dw e− 1 2G 2w2+Hw qθ(w) (5)
Notice that for common deterministic activations such as ReLU and sign, the function ϕ has analytic and smooth expressions (see Appendix A.8). The same holds for the function ψ when qθ(w) is Gaussian (continuous weights) or a mixture of atoms (discrete weights). At the last layer we impose PL+1(y|z) = I(y = sign(z)) in binary classification tasks and PL+1(y|z) = I(y = arg max(z)) in multi-class classification (see Appendix A.9). While in our experiments we use hard constraints for the final output, therefore solving a constraint satisfaction problem, it would be interesting to also consider soft constraints and introduce a temperature, but this is beyond the scope of our work.
Start and end of message passing. At the beginning of a new PasP iteration t, we reset the messages (see Appendix A) and run message passing for τmax iterations. We then compute the new prior’s parameters θt+1 from the posterior given by the message passing.
BP Forward pass. After initialization of the messages at time τ = 0, for each following time we propagate a set of message from the first to the last layer and then another set from the last to the first. For an intermediate layer ` the forward pass reads
x̂`,τin→k = ∂Bϕ ` ( B`,τ−1in→k , A `,τ−1 in , ω `−1,τ in , V `−1,τ in ) (6)
∆`,τin = ∂ 2 Bϕ ` ( B`,τ−1in , A `,τ−1 in , ω `−1,τ in , V `−1,τ in ) (7)
m`,τki→n = ∂Hψ(H `,τ−1 ki→n , G `,τ−1 ki , θ ` ki) (8)
σ`,τki = ∂ 2 Hψ(H `,τ−1 ki , G `,τ−1 ki , θ ` ki) (9)
V `,τkn = ∑ i (( m`,τki→n )2 ∆`,τin + σ `,τ ki (x̂ `,τ in→k) 2 + σ`,τki ∆ `,τ in ) (10)
ω`,τkn→i = ∑ i′ 6=i m`,τki′→nx̂ `,τ i′n→k (11)
The equations for the first layer differ slightly and in an intuitive way from the ones above (see Appendix A.3).
BP Backward pass. The backward pass updates a set of messages from the last to the first layer:
g`,τkn→i = ∂ωϕ `+1 ( B`+1,τkn , A `+1,τ kn , ω `,τ kn→i, V `,τ kn ) (12)
Γ`,τkn = −∂ 2 ωϕ `+1 ( B`+1,τkn , A `+1,τ kn , ω `,τ kn , V `,τ kn ) (13)
A`,τin = ∑ k ( (m`,τki→n) 2 + σ`,τki ) Γ`,τkn − σ `,τ ki ( g`,τkn→i )2 (14)
B`,τin→k = ∑ k′ 6=k m`,τk′i→ng `,τ k′n→i (15)
G`,τki = ∑ n ( (x̂`,τin→k) 2 + ∆`,τin ) Γ`,τkn −∆ `,τ in ( g`,τkn→i )2 (16)
H`,τki→n = ∑ n′ 6=n x̂`,τin′→kg `,τ kn′→i (17)
As with the forward pass, we add the caveat that for the last layer the equations are slightly different from the ones above.
Computational complexity The message passing equations boil down to element-wise operations and tensor contractions that we easily implement using the GPU friendly julia library Tullio.jl (Abbott et al., 2021). For a layer of input and output size N and considering a batch-size of B, the time complexity of a forth-and-back iteration is O(N2B) for all message passing algorithms (BP, BPI, MF, and AMP), the same as SGD. The prefactor varies and it is generally larger than SGD (see Appendix B.9). Also, time complexity for message passing is proportional to τmax (which we typically set to 1). We provide our implementation in the GitHub repo anonymous.
Algorithm 1: BP for deep neural networks // Message passing used in the PasP Eq. 3 to approximate. // the mini-batch posterior. // Here we specifically refer to BP updates. // BPI, MF, and AMP updates take the same form but using // the rules in Appendix A.4, A.5, and A.7 respectively
1 Initialize messages. 2 for τ = 1, . . . τmax do // Forward Pass 3 for l = 0, . . . , L do 4 compute x̂`,∆` using (6, 7) 5 computem`,σ` using (8, 9) 6 compute V`,ω` using (10, 11)
// Backward Pass 7 for l = L, . . . , 0 do 8 compute g`,Γ` using (12, 13) 9 computeA`,B` using (14, 15)
10 computeG`,H` using (16, 17)

heading: 4 NUMERICAL RESULTS
text: We implement our message passing algorithms on neural networks with continuous and binary weights and with binary activations. In our experiments we fix τmax = 1. We typically do not observe an increase in performance taking more steps, except for some specific cases and in particular for MF layers. We remark that for τmax = 1 the BP and the BPI equations are identical, so in most of the subsequent numerical results we will only investigate BP.
We compare our algorithms with a SGD-based algorithm adapted to binary architectures (Hubara et al., 2016) which we call BinaryNet along the paper (see Appendix B.6 for details). Comparison of Bayesian predictions are with the gradient-based Expectation Backpropagation (EBP) algorithm (Soudry et al., 2014a), also able to deal with discrete weights and activations. In all architectures we avoid the use of bias terms and batch-normalization layers.
We find that message-passing algorithms are able to train generic MLP architectures with varying numbers and sizes of hidden layers. As for the datasets, we are able to perform both binary classification and multi-class classification on standard computer vision datasets such as MNIST, Fashion-MNIST, and CIFAR-10. Since these datasets consist of 10 classes, for the binary classification task we divide each dataset in two classes (even vs odd).
We report that message passing algorithms are able to solve these optimization problems with generalization performance comparable to or better than SGD-based algorithms. Some of the message passing algorithms (BP and AMP in particular) need fewer epochs to achieve low error than the ones required by SGD-based algorithms, even if adaptive methods like Adam are considered. Timings of our GPU implementations of message passing algorithms are competitive with SGD (see Appendix B.9).

heading: 4.1 EXPERIMENTS ACROSS ARCHITECTURES
text: We select a specific task, multi-class classification on Fashion-MNIST, and we compare the message passing algorithms with BinaryNet for different choices of the architecture (i.e. we vary the number and the size of the hidden layers). In Fig.1 (Left) we present the learning curves for a MLP with 3 hidden layers with 501 units with binary weights and activations. Similar results hold in our experiments with 2 or 3 hidden layers of 101, 501 or 1001 units and with batch sizes from 1 to from 1024. The parameters used in our simulations are reported in Appendix B.3. Results on networks with continuous weights can be found in Fig.2 (Right).

heading: 4.2 SPARSE LAYERS
text: Since the BP algorithm has notoriously been successful on sparse graphs, we perform a straightforward implementation of pruning at initialization, i.e. we impose a random boolean mask on the weights that we keep fixed along the training. We call sparsity the fraction of zeroed weights. This kind of non-adaptive pruning is known to largely hinder learning (Frankle et al., 2021; Sung et al., 2021). In the right panel of Fig. 1, we report results on sparse binary networks in which we train a MLP with 2 hidden layers of 101 units on the MNIST dataset. For reference, results on pruning quantized/binary networks can be found in Refs. (Han et al., 2016; Ardakani et al., 2017; Tung & Mori, 2018; Diffenderfer & Kailkhura, 2021). Experimenting with sparsity up to 90%, we observe that BP and MF perform better than BinaryNet. AMP struggles behind BinaryNet instead.

heading: 4.3 EXPERIMENTS ACROSS DATASETS
text: We now fix the architecture, a MLP with 2 hidden layers of 501 neurons each with binary weights and activations. We vary the dataset, i.e. we test the BP-based algorithms on standard computer vision benchmark datasets such as MNIST, Fashion-MNIST and CIFAR-10, in both the multi-class and binary classification tasks. In Tab. 1 we report the final test errors obtained by the message passing algorithms compared to the BinaryNet baseline. See Appendix B.4 for the corresponding training errors and the parameters used in the simulations. We mention that while the test performance is mostly comparable, the train error tends to be lower for the message passing algorithms.

heading: 4.4 LOCALLY BAYESIAN ERROR
text: The message passing framework used as an estimator of the mini-batch posterior marginals allows us to perform approximate Bayesian prediction, i.e. averaging the pointwise predictions over the approximate posterior. We observe better generalization error from Bayesian predictions compared to point-wise ones, showing that the marginals retain useful information. However, we roughly estimate the marginals with the PasP mini-batch procedure (the exact ones should be computed with a full-batch procedure, but this converges with difficulty in our tests). Since BP-based algorithms tend to focus on dense states (as also confirmed by the local energy measure performed in Appendix B.5), the Bayesian error we compute can be considered as a local approximation of the full one. We report results for binary classification on the MNIST dataset in Fig. 2, and we observe the same performance increase on different datasets and architectures. We obtain the Bayesian prediction from the output marginal given by a single forward pass of the message passing. To obtain good Bayesian estimates it is important that the posterior distribution does not concentrate too much, otherwise the Bayesian prediction will converge to the prediction of a single configuration.
In Fig.2 we also perform a comparison of BP (point-wise and Bayesian) with SGD and another algorithm able to perform Bayesian predictions, Expectation Backpropagation (Soudry et al., 2014a) see Appendix B.7 for implementation details.

heading: 4.5 CONTINUAL LEARNING
text: Given the high local entropy (i.e. the flatness) of the solutions found by the BP-based algorithms (see Appendix B.5), we perform additional tests in a classic setting, continual learning, where the
possibility of locally rearranging the solutions while keeping low training error can be an advantage. When a deep network is trained sequentially on different tasks, it tends to forget exponentially fast previously seen tasks while learning new ones (McCloskey & Cohen, 1989; Robins, 1995; Fusi et al., 2005). Recent work (Feng & Tu, 2021) has shown that searching for a flat region in the loss landscape can indeed help to prevent catastrophic forgetting. Several heuristics have been proposed to mitigate the problem (Kirkpatrick et al., 2017; Aljundi et al., 2018; Zenke et al., 2017; Laborieux et al., 2021) but all require specialized adjustments to the loss or the dynamics .
Here we show instead that our message passing schemes are naturally prone to learn multiple tasks sequentially, mitigating the characteristic memory issues of the gradient-based schemes without the need for explicit modifications. As a prototypical experiment, we sequentially trained a multi-layer neural network on 6 different versions of the MNIST dataset, where the pixels of the images have been randomly permuted (Goodfellow et al., 2013), giving a fixed budget of 40 epochs on each task. We present the results for a two hidden layer neural network with 2001 units on each layer (see Appendix B.3 for details). As can be seen in Fig. 3, at the end of the training the BP algorithm is able to reach good generalization performances on all the tasks. We compared the BP performance with BinaryNet, which already performs better than SGD with continuous weights (see the discussion in Laborieux et al. (2021)). While our BP implementation is not competitive with ad-hoc techniques specifically designed for this problem, it beats non-specialized heuristics. Moreover, we believe that specialized approaches like the one of Laborieux et al. (2021) can be adapted to message passing as well.

heading: 5 DISCUSSION AND CONCLUSIONS
text: While successful in many fields, message passing algorithms, have notoriously struggled to scale to deep neural networks training problems. Here we have developed a class of fBP-based message passing algorithms and used them within an update scheme, Posterior-as-Prior (PasP), that makes it possible to train deep and wide multilayer perceptrons by message passing.
We performed experiments binary activations and either binary or continuous weights. Future work should try to include different activations, biases, batch-normalization, and convolutional layers as well. Another interesting direction is the algorithmic computation of the (local) entropy of the model from the messages.
Further theoretical work is needed for a more complete understanding of the robustness of our methods. Recent developments in message passing algorithms (Rangan et al., 2019) and related theoretical analysis (Goldt et al., 2020) could provide fruitful inspirations. While our algorithms can be used for approximate Bayesian inference, exact posterior calculation is still out of reach for message passing approaches and much technical work is needed in that direction.

heading: Appendices

heading: CONTENTS

heading: A BP-based message passing algorithms 14
text: A.1 Preliminary considerations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
A.2 Derivation of the BP equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
A.3 BP equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
A.4 BPI equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
A.5 MF equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
A.6 Derivation of the AMP equations . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
A.7 AMP equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
A.8 Activation Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
A.9 The ArgMax layer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

heading: B Experimental details 25
text: B.1 Hyper-parameters of the BP-based scheme . . . . . . . . . . . . . . . . . . . . . . 25
B.2 Damping scheme for the message passing . . . . . . . . . . . . . . . . . . . . . . 25
B.3 Architectures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
B.4 Varying the dataset . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
B.5 Local energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
B.6 SGD implementation (BinaryNet) . . . . . . . . . . . . . . . . . . . . . . . . . . 28
B.7 EBP implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
B.8 Unit polarization and overlaps . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
B.9 Computational performance: varying batch-size . . . . . . . . . . . . . . . . . . . 30

heading: A BP-BASED MESSAGE PASSING ALGORITHMS

heading: A.1 PRELIMINARY CONSIDERATIONS
text: Given a mini-batch B = {(xn, yn)}n, the factor graph defined by Eqs. (1, 2, 18) is explicitly written as:
P (W,x1:L | B, θ) ∝ L∏ `=0 ∏ k,n P `+1 ( x`+1kn ∣∣∣∣ ∑ i W `kix ` in ) ∏ k,i,` qθ(W ` ki), (18)
where x0n = xn, x L+1 n = yn. The derivation of the BP equations for this model is straightforward albeit lengthy and involved. It is obtained following the steps presented in multiple papers, books, and reviews, see for instance (Mézard & Montanari, 2009; Zdeborová & Krzakala, 2016; Mézard, 2017), although it has not been attempted before in deep neural networks. It should be noted that a (common) approximation that we take here with respect to the standard BP scheme, is that messages are assumed to be Gaussian distributed and therefore parameterized by their mean and variance. This goes by the name of relaxed belied propagation (rBP), just referred to as BP throughout the paper.
We derive the BP equations in A.2 and present them all together in A.3. From BP, we derive other 3 message passing algorithms useful for the deep network training setting, all of which are well known to the literature: BP-Inspired (BPI) message passing A.4, mean-field (MF) A.5, and approximate
message passing (AMP) A.7. The AMP derivation is the more involved and given in A.6. In all these cases, message updates can be divided in a forward pass and a backward pass, as also done in Fletcher et al. (2018) in a multi-layer inference setting. The BP algorithm is compactly reported in Algorithm 1.
In our notation, ` denotes the layer index, τ the BP iteration index, k an output neuron index, i an input neuron index, and n a sample index.
We report below, for convenience, some of the considerations also present in the main text.
Meaning of messages. All the messages involved in the message passing equations can be understood in terms of cavity marginals or full marginals (as mentioned in the introduction BP is also known as the Cavity Method, see Mézard & Montanari (2009)). Of particular relevance are the quantities m`ki and σ ` ki, denoting the mean and variance of the weights W ` ki. The quantities x̂ ` in and ∆`in instead denote mean and variance of the i-th neuron’s activation in layer ` in correspondence of an input xn.
Scalar free energies. All message passing schemes can be expressed using the following scalar functions, corresponding to single neuron and single weight effective free-energies respectively:
ϕ`(B,A, ω, V ) = log ∫ dx dz e− 1 2Ax 2+Bx P ` (x | z) e− (ω−z)2 2V , (19)
ψ(H,G, θ) = log ∫ dw e− 1 2G 2w2+Hw qθ(w). (20)
These free energies will naturally arise in the derivation of the BP equations in Appendix A.2. For the last layer, the neuron function has to be slightly modified:
ϕL+1(y, ω, V ) = log ∫ dz PL+1 (y | z) e− (ω−z)2 2V . (21)
Notice that for common deterministic activations such as ReLU and sign, the function ϕ has analytic and smooth expressions that we give in Appendix A.8. Same goes for ψ when qθ(w) is Gaussian (continuous weights) or a mixture of atoms (discrete weights). At the last layer we impose PL+1(y|z) = I(y = sign(z)) in binary classification tasks. For multi-class classification instead, we have to adapt the formalism to vectorial pre-activations z and assume PL+1(y|z) = I(y = arg max(z)) (see Appendix A.9). While in our experiments we use hard constraints for the final output, therefore solving a constraint satisfaction problem, it would be interesting to also consider generic loss functions. That would require minimal changes to our formalism, but this is beyond the scope of our work.
Binary weights. In our experiments we use ±1 weights in each layer. Therefore each marginal can be parameterized by a single number and our prior/posterior takes the form
qθ(W ` ki) ∝ eθ
` kiW ` ki (22)
The effective free energy function Eq. 20 becomes
ψ(H,G, θ`ki) = log 2 cosh(H + θ ` ki) (23)
and the messages G can be dropped from the message passing.
Start and end of message passing. At the beginning of a new PasP iteration t, we reset the messages to zero and run message passing for τmax iterations. We then compute the new prior qθt+1(W) from the posterior given by the message passing iterations.

heading: A.2 DERIVATION OF THE BP EQUATIONS
text: In order to derive the BP equations, we start with the following portion of the factor graph reported in Eq. 18 in the main text, describing the contribution of a single data example in the inner loop of the PasP updates:
L∏ `=0 ∏ k P `+1 ( x`+1kn ∣∣∣∣ ∑ i W `kix ` in ) where x0n = xn, x L+1 n = yn. (24)
where we recall that the quantity x`kn corresponds to the activation of neuron k in layer ` in correspondence of the input example n.
Let us start by analyzing the single factor:
P `+1 ( x`+1kn ∣∣∣∣ ∑ i W `kix ` in ) (25)
We refer to messages that travel from input to output in the factor graph as upgoing or upwards messages, while to the ones that travel from output to input as downgoing or backwards messages.
Factor-to-variable-W messages The factor-to-variable-W messages read:
ν̂`+1kn→ki(W ` ki) ∝ ∫ ∏ i′ 6=i dν`ki′→n(W ` ki′) ∏ i′ dν`i′n→k(x ` i′n) dν↓(x `+1 kn ) P `+1 ( x`+1kn ∣∣∣∣ ∑ i′ W `ki′x ` i′n ) (26)
where ν↓ denotes the messages travelling downwards (from output to input) in the factor graph.
We denote the means and variances of the incoming messages respectively with m`ki→n, x̂ ` in→k and σ`ki→n, ∆ ` in→k:
m`ki→n = ∫ dν`ki→n(W ` ki) W ` ki (27)
σ`ki→n = ∫ dν`ki→n(W ` ki) ( W `ki −m`ki→n )2 (28)
x̂`in→k = ∫ dν`in→k(x ` in) x ` in (29)
∆`in→k = ∫ dν`in→k(x ` in) ( x`in − x̂`in→k )2 (30)
We now use the central limit theorem to observe that with respect to the incoming messages distributions - assuming independence of these messages - in the large input limit the preactivation is a Gaussian random variable: ∑
i′ 6=i
W `ki′x ` i′n ∼ N (ω`kn→i, V `kn→i) (31)
where:
ω`kn→i = Eν ∑ i′ 6=i W `ki′x ` i′n  = ∑ i′ 6=i m`ki′→n x̂ ` i′n→k (32)
V `kn→i = V arν ∑ i′ 6=i W `ki′x ` i′n  = ∑ i′ 6=i ( σ`ki′→n ∆ ` i′n→k + ( m`ki′→n )2 ∆`i′n→k + σ ` ki′→n ( x̂`i′n→k )2) (33)
Therefore we can rewrite the outgoing messages as:
ν̂`+1kn→i(W ` ki) ∝ ∫ dz dν`in→k(x ` in) dν↓(x `+1 kn ) e − ( z−ωkn→i−W ` kix ` in) 2 2Vkn→i P `+1 ( x`+1kn ∣∣∣∣ z) (34) We now assume W `kix ` in to be small compared to the other terms. With a second order Taylor expansion we obtain:
ν̂`kn→i(W ` ki) ∝ ∫ dz dν↓(x `+1 kn ) e − (z−ωkn→i) 2 2Vkn→i P `+1 ( x`+1kn ∣∣∣∣ z) × ( 1 +
z − ωkn→i Vkn→i x̂`in→kW ` ki + (z − ωkn→i)2 − Vkn→i 2Vkn→i
( ∆ + ( x̂`in→k )2) ( W `ki )2) (35)
Introducing now the function:
ϕ`(B,A, ω, V ) = log ∫ dx dz e− 1 2Ax 2+Bx P ` (x|z) e− (ω−z)2 2V (36)
and defining:
g`kn→i = ∂ωϕ `+1(B`+1, A`+1, ω`kn→i, V ` kn→i) (37)
Γ`kn→i = −∂2ωϕ`+1(B`+1, A`+1, ω`kn→i, V `kn→i) (38)
the expansion for the log-message reads:
log ν̂`kn→i(W ` ki) ≈ const+ x̂`in→k g`kn→iW `ki
− 1 2
(( ∆`in→k + ( x̂`in→k )2) Γ`kn→i −∆`in→k ( g`kn→i )2) ( W `ki )2 (39)
Factor-to-variable-x messages The derivation of these messages is analogous to the factor-tovariable-W ones in Eq. 26 just reported. The final result for the log-message is:
log ν̂`kn→i(x ` in) ≈const+m`ki→n g`kn→ix`in
− 1 2
(( σ`ki→n + ( m`ki→n )2) Γ`kn→i − σ`ki→n ( g`kn→i )2) ( x`in )2
(40)
Variable-W-to-output-factor messages The message from variable W `ki to the output factor kn reads:
ν`ki→n(W ` ki) ∝ P `θki(W ` ki)e
∑ n′ 6=n log ν̂ ` kn′→i(W ` ki)
≈ P `θki(W ` ki)e
H`ki→nW ` ki− 12G ` ki→n(W ` ki)
2
(41)
where we have defined:
H`ki→n = ∑ n′ 6=n x̂`in′→k g ` kn′→i (42)
G`ki→n = ∑ n′ 6=n (( ∆`in′→k + ( x̂`in′→k )2) Γ`kn′→i −∆`in′→k ( g`kn′→i )2) (43)
Introducing now the effective free energy:
ψ(H,G, θ) = log ∫ dW P `θ (W ) e HW− 12GW 2
(44)
we can express the first two cumulants of the message ν`ki→n(W ` ki) as:
m`ki→n = ∂Hψ(H ` ki→n, G ` ki→n, θki) (45)
σ`ki→n = ∂ 2 Hψ(H ` ki→n, G ` ki→n, θki) (46)
Variable-x-to-input-factor messages We can write the downgoing message as:
ν↓(x ` in) ∝ e
∑ k log ν̂ ` kn→i(x ` in)
≈ eB ` inx− 12A ` inx 2
(47)
where:
B`in = ∑ n m`ki→n g ` kn→i (48)
A`in = ∑ n (( σ`ki→n + ( m`ki→n )2) Γ`kn→i − σ`ki→n ( g`+1kn→i )2) (49)
Variable-x-to-output-factor messages By defining the following cavity quantities:
B`in→k = B ` in→k −m`ki→n g`kn→i (50)
A`in→k = A ` in→k − (( σ`ki→n + ( m`ki→n )2) Γ`kn→i − σ`ki→n ( g`kn→i )2) (51)
and the following non-cavity ones:
ω`kn = ∑ i m`ki→n x̂ ` in→k (52)
V `kn = ∑ i ( σ`ki→n ∆ ` in→k + ( m`ki→n )2 ∆`in→k + σ ` ki→n ( x̂`i′n→k )2) (53)
we can express the first 2 cumulants of the upgoing messages as:
x̂`in→k = ∂Bϕ `(B`in→k, A ` in→k, ω `−1 in , V `−1 in ) (54)
∆`in→k = ∂ 2 Bϕ `(B`in→k, A ` in→k, ω `−1 in , V `−1 in ) (55)
Wrapping it up Additional but straightforward considerations are required for the final input and output layers (` = 0 and ` = L respectively), since they do not receive messages from below and above respectively. In the end, thanks to independence assumptions and the central limit theorem that we used throughout the derivations, we arrive to a closed set of equations involving the means and the variances (or otherwise the corresponding natural parameters) of the messages. Within the same approximation assumption, we also replace the cavity quantities corresponding to variances with the non-cavity counterparts. Dividing the update equations in a forward and backward pass, and ordering them using time indexes in such a way that we have an efficient flow of information, we obtain the set of BP equations presented in the main text Eqs. (6-17) and in the Appendix Eqs. (60-71).

heading: A.3 BP EQUATIONS
text: We report here the end result of the derivation in last section, the complete set of BP equations also presented in the main text as Eqs. (6-17).
Initialization At τ = 0:
B`,0in→k = 0 (56)
A`,0in = 0 (57)
H`,0ki→n = 0 (58)
G`,0ki = 0 (59)
Forward Pass At each τ = 1, . . . , τmax, for ` = 0, . . . , L:
x̂`,τin→k = ∂Bϕ `(B`,τ−1in→k , A `,τ−1 in , ω `−1,τ in , V `−1,τ in ) (60)
∆`,τin = ∂ 2 Bϕ `(B`,τ−1in , A `,τ−1 in , ω `−1,τ in , V `−1,τ in ) (61)
m`,τki→n = ∂Hψ(H `,τ−1 ki→n , G `,τ−1 ki , θ ` ki) (62)
σ`,τki = ∂ 2 Hψ(H `,τ−1 ki , G `,τ−1 ki , θ ` ki) (63)
V `,τkn = ∑ i (( m`,τki )2 ∆`,τin + σ `,τ−1 ki ( x̂`,τi′n )2 + σ`,τ−1ki ∆ `,τ in ) (64)
ω`,τkn→i = ∑ i′ 6=i m`,τki′→n x̂ `,τ i′n→k (65)
In these equations for simplicity we abused the notation, in fact for the first layer x̂`=0,τn is fixed and given by the input xn while ∆`=0,τn = 0 instead.
Backward Pass For τ = 1, . . . , τmax, for ` = L, . . . , 0 :
g`,τkn→i = ∂ωϕ `+1(B`+1,τkn , A `+1,τ kn , ω `,τ kn→i, V `,τ kn ) (66)
Γ`,τkn = −∂ 2 ωϕ `+1(B`+1,τkn , A `+1,τ kn , ω `,τ kn , V `,τ kn ) (67) A`,τin = ∑ k ((( m`,τki )2 + σ`,τki ) Γ`,τkn − σ `,τ ki ( g`,τkn )2) (68)
B`,τin→k = ∑ k′ 6=k m`,τk′i→n g `,τ k′n→i (69)
G`,τki = ∑ n ((( x̂`,τin )2 + ∆`,τin ) Γ`,τkn −∆ `,τ in ( g`,τkn )2) (70)
H`,τki→n = ∑ n′ 6=n x̂`,τin′→k g `,τ kn′→i (71)
In these equations as well we abused the notation: calling L the number of hidden neuron layers, when ` = L one should use ϕL+1(y, ω, V ) from Eq. 21 instead of ϕL+1(B,A, ω, V ).

heading: A.4 BPI EQUATIONS
text: The BP-Inspired algorithm (BPI) is obtained as an approximation of BP replacing some cavity quantities with their non-cavity counterparts. What we obtain is a generalization of the single layer algorithm of Baldassi et al. (2007).
Forward pass.
x̂`,τin = ∂Bϕ ` ( B`,τ−1in , A `,τ−1 in , ω `−1,τ in , V `−1,τ in ) (72)
∆`,τin = ∂ 2 Bϕ ` ( B`,τ−1in , A `,τ−1 in , ω `−1,τ in , V `−1,τ in ) (73)
m`,τki = ∂Hψ(H `,τ−1 ki , G `,τ−1 ki , θ ` ki) (74)
σ`,τki = ∂ 2 Hψ(H `,τ−1 ki , G `,τ−1 ki , θ ` ki) (75)
V `,τkn = ∑ i (( m`,τki )2 ∆`,τin + σ `,τ ki (x̂ `,τ in ) 2 + σ`,τki ∆ `,τ in ) (76)
ω`,τkn = ∑ i m`,τki x̂ `,τ in (77)
Backward pass.
g`,τkn→i = ∂ωϕ `+1 ( B`+1,τkn , A `+1,τ kn , ω `,τ kn −m `,τ ki x̂ `,τ ai , V `,τ kn ) (78)
Γ`,τkn = −∂ 2 ωϕ `+1 ( B`+1,τkn , A `+1,τ kn , ω `,τ kn , V `,τ kn ) (79)
A`,τin = ∑ k ( (m`,τki ) 2 + σ`,τki ) Γ`,τkn − σ `,τ ki ( g`,τkn )2 (80)
B`,τin = ∑ k m`,τki g `,τ kn→i (81)
G`,τki = ∑ n ( (x̂`,τin ) 2 + ∆`,τin ) Γ`,τkn −∆ `,τ in ( g`,τkn )2 (82)
H`,τki = ∑ n x̂`,τin g `,τ kn→i (83)

heading: A.5 MF EQUATIONS
text: The mean-field (MF) equations are obtained as a further simplification of BPI, using only non-cavity quantities. Although the simplification appears minimal at this point, we empirically observe a non-negligible discrepancy between the two algorithms in terms of generalization performance and computational time.
Forward pass.
x̂`,τin = ∂Bϕ ` ( B`,τ−1in , A `,τ−1 in , ω `−1,τ in , V `−1,τ in ) (84)
∆`,τin = ∂ 2 Bϕ ` ( B`,τ−1in , A `,τ−1 in , ω `−1,τ in , V `−1,τ in ) (85)
m`,τki = ∂Hψ(H `,τ−1 ki , G `,τ−1 ki , θ ` ki) (86)
σ`,τki = ∂ 2 Hψ(H `,τ−1 ki , G `,τ−1 ki , θ ` ki) (87)
V `,τkn = ∑ i (( m`,τki )2 ∆`,τin + σ `,τ ki (x̂ `,τ in ) 2 + σ`,τki ∆ `,τ in ) (88)
ω`,τkn = ∑ i m`,τki x̂ `,τ in (89)
Backward pass.
g`,τkn = ∂ωϕ `+1 ( B`+1,τkn , A `+1,τ kn , ω `,τ kn , V `,τ kn ) (90)
Γ`,τkn = −∂ 2 ωϕ `+1 ( B`+1,τkn , A `+1,τ kn , ω `,τ kn , V `,τ kn ) (91)
A`,τin = ∑ k ( (m`,τki ) 2 + σ`,τki ) Γ`,τkn − σ `,τ ki ( g`,τkn )2 (92)
B`,τin = ∑ k m`,τki g `,τ kn (93)
G`,τki = ∑ n ( (x̂`,τin ) 2 + ∆`,τin ) Γ`,τkn −∆ `,τ in ( g`,τkn )2 (94)
H`,τki = ∑ n x̂`,τin g `,τ kn (95)

heading: A.6 DERIVATION OF THE AMP EQUATIONS
text: In order to obtain the AMP equations, we approximate cavity quantities with non-cavity ones in the BP equations Eqs. (60-71) using a first order expansion. We start with the mean activation:
x̂`,τin→k =∂Bϕ `(B`,τ−1in −m `,τ−1 ki→n g `,τ−1 kn→i , A `,τ−1 in , ω `−1,τ in , V `−1,τ in )
≈∂Bϕ`(B`,τ−1in , A `,τ−1 in , ω `−1,τ in , V `−1,τ in )
−m`,τ−1ki→n g `,τ−1 kn→i∂ 2 Bϕ `(B`,τ−1in , A `,τ−1 in , ω `−1,τ in , V `−1,τ in )
≈x̂`,τin −m `,τ−1 ki g `,τ−1 kn ∆ `,τ in (96)
Analogously, for the weight’s mean we have:
m`,τki→n = ∂Hψ(H `,τ−1 ki − x̂ `,τ−1 in→k g `,τ−1 kn→i , G `,τ−1 ki , θ ` ki)
≈ ∂Hψ(H`,τ−1ki , G `,τ−1 ki , θ ` ki)− x̂ `,τ−1 in→k g `,τ−1 kn→i∂ 2 Hψ(H `,τ−1 ki , G `,τ−1 ki , θ ` ki) ≈ m`,τki − x̂ `,τ−1 in g `,τ−1 kn σ `,τ ki . (97)
This brings us to:
ω`,τkn = ∑ i m`,τki→n x̂ `,τ in→k
≈ ∑ i m`,τki x̂ `,τ in − g `,τ−1 kn ∑ i σ`,τki x̂ `,τ in x̂ `,τ−1 in − g `,τ−1 kn ∑ i m`,τki m `,τ−1 ki ∆ `,τ in
+ (g`,τ−1kn ) 2 ∑ i σ`,τki m `,τ−1 ki x̂ `,τ−1 in ∆ `,τ in (98)
Let us now apply the same procedure to the other set of cavity messages:
g`,τkn→i =∂ωϕ `+1(B`+1,τkn , A `+1,τ kn , ω `,τ kn −m `,τ ki→n x̂ `,τ in→k, V `,τ kn )
≈∂ωϕ`+1(B`+1,τkn , A `+1,τ kn , ω `,τ kn , V `,τ kn )
−m`,τki→n x̂ `,τ in→k∂ 2 ωϕ `+1(B`+1,τkn , A `+1,τ kn , ω `,τ kn , V `,τ kn )
≈g`,τkn +m `,τ ki x̂ `,τ in Γ `,τ kn (99)
B`,τin = ∑ k m`,τki→n g `,τ kn→i
≈ ∑ k m`,τki g `,τ kn − x̂in ∑ k ( g`,τkn )2 σ`,τki + x̂ `,τ in ∑ k (m`,τki ) 2Γ`,τkn
− (x̂`,τin ) 2 ∑ k σ`,τki m `,τ ki g `,τ knΓ `,τ kn (100)
H`,τki = ∑ n x̂`,τin→k g `,τ kn→i
≈ ∑ n x̂`,τin g `,τ kn +m `,τ ki ∑ n ( x̂`,τin )2 Γ`,τkn −m `,τ ki ∑ n (g`,τkn ) 2∆`,τin
− (m`,τki ) 2 ∑ n g`,τknΓ `,τ kn∆ `,τ in x̂ `,τ in (101)
We are now able to write down the full AMP equations, that we present in the next section.

heading: A.7 AMP EQUATIONS
text: In summary, in the last section we derived the AMP algorithm as a closure of the BP messages passing over non-cavity quantities, relying on some statistical assumptions on messages and interactions. With respect to the MF message passing, we find some additional terms that go under the name of Onsager corrections. In-depth overviews of the AMP (also known as Thouless-Anderson-Palmer (TAP)) approach can be found in Refs. (Zdeborová & Krzakala, 2016; Mézard, 2017; Gabrié, 2020). The final form of the AMP equations for the multi-layer perceptron is given below.

heading: Initialization At τ = 0:
text: B`,0in = 0 (102)
A`,0in = 0 (103)
H`,0ki = 0 or some values (104)
G`,0ki = 0 or some values (105)
g`,0kn = 0 (106)
Forward Pass At each τ = 1, . . . , τmax, for ` = 0, . . . , L:
x̂`,τin =∂Bϕ `(B`,τ−1in , A `,τ−1 in , ω `−1,τ in , V `−1,τ in ) (107)
∆`,τin =∂ 2 Bϕ `(B`,τ−1in , A `,τ−1 in , ω `−1,τ in , V `−1,τ in ) (108) m`,τki =∂Hψ(H `,τ−1 ki , G `,τ−1 ki , θ ` ki) (109)
σ`,τki =∂ 2 Hψ(H `,τ−1 ki , G `,τ−1 ki , θ ` ki) (110)
V `,τkn = ∑ i (( m`,τki )2 ∆`,τin + σ `,τ ki ( x̂`,τi′n )2 + σ`,τki ∆ `,τ in ) (111)
ω`,τkn = ∑ i m`,τki x̂ `,τ in − g `,τ−1 kn ∑ i σ`,τki x̂ `,τ in x̂ `,τ−1 in − g `,τ−1 kn ∑ i m`,τki m `,τ−1 ki ∆ `,τ in
+ (g`,τ−1kn ) 2 ∑ i σ`,τki m `,τ−1 ki x̂ `,τ−1 in ∆ `,τ in (112)
Backward Pass
g`,τkn =∂ωϕ `+1(B`+1,τkn , A `+1,τ kn , ω `,τ kn→i, V `,τ kn ) (113)
Γ`,τkn =− ∂ 2 ωϕ `+1(B`+1,τkn , A `+1,τ kn , ω `,τ kn , V `,τ kn ) (114) A`,τin = ∑ k ((( m`,τki )2 + σ`,τki ) Γ`,τkn − σ `,τ ki ( g`,τkn )2) (115)
B`,τin = ∑ k m`,τki g `,τ kn − x̂in ∑ k ( g`,τkn )2 σ`,τki + x̂ `,τ in ∑ k (m`,τki ) 2Γ`,τkn
− (x̂`,τin ) 2 ∑ k σ`,τki m `,τ ki g `,τ knΓ `,τ kn (116)
G`,τki = ∑ n ((( x̂`,τin )2 + ∆`,τin ) Γ`,τkn −∆ `,τ in ( g`,τkn )2) (117)
H`,τki = ∑ n x̂`,τin g `,τ kn +m `,τ ki ∑ n ( x̂`,τin )2 Γ`,τkn −m `,τ ki ∑ n (g`,τkn ) 2∆`,τin
− (m`,τki ) 2 ∑ n g`,τknΓ `,τ kn∆ `,τ in x̂ `,τ in (118)

heading: A.8 ACTIVATION FUNCTIONS
text: A.8.1 SIGN
In most of our experiments we use sign activations in each layer. With this choice, the neuron’s free energy 19 takes the form
ϕ(B,A, ω, V ) = log 1 2 ∑ x∈{−1,+1} eBxH ( − xω√ V )+ 1 2 log(2πV ), (119)
where
H = 1 2 erfc ( x√ 2 ) . (120)
Notice that for sign activations the messages A can be dropped.
A.8.2 RELU
For ReLU(x) = max(0, x) activations the free energy 19 becomes
ϕ(B,A, ω, V ) = ∫ dxdz e− 1 2Ax 2+Bx δ(x−max(0, z)) e− (ω−z)2 2V (121)
= log ( H ( ω√ V ) + N (ω;B/A, V + 1A ) AN (B; 0, A) H ( − BV + ω√ V +AV 2 )) + 1 2 log(2πV ),
(122)
where
N (x;µ,Σ) = 1√ 2πΣ e− (x−µ)2 2Σ . (123)

heading: A.9 THE ARGMAX LAYER
text: In order to perform multi-class classification, we have to perform an argmax operation on the last layer of the neural network. Call zk, for k = 1, . . . ,K, the Gaussian random variables output of the last layer of the network in correspondence of some input x. Assuming the correct label is class k∗, the effective partition function Zk∗ corresponding to the output constraint reads:
Zk∗ = ∫ ∏ k dzkN (zk;ωk, Vk) ∏ k 6=k∗ Θ(zk∗ − zk) (124)
= ∫ dzk∗ N (zk∗ ;ωk∗ , Vk∗) ∏ k 6=k∗ H ( −zk ∗ − ωk√ Vk ) (125)
Here Θ(x) is the Heaviside indicator function and we used the definition of H from Eq. 120. The integral on the last line cannot be expressed analytically, therefore we have to resort to approximations.

heading: A.9.1 APPROACH 1: JENSEN INEQUALITY
text: Using the Jensen inequality we obtain:
φk∗ = logZk∗ = logEz∼N (ωk∗ ,Vk∗ ) ∏ k 6=k∗ H ( −z − ωk√ Vk ) (126)
≥ ∑ k 6=k∗ Ez∼N (ωk∗ ,Vk∗ ) logH ( −z − ωk√ Vk ) (127)
Reparameterizing the expectation we have:
φ̃k∗ = ∑ k 6=k∗ E ∼N (0,1) logH ( −ωk ∗ + √ Vk∗ − ωk√ Vk ) (128)
The derivative ∂ωk φ̃k∗ and ∂ 2 ωk φ̃k∗ that we need can then be estimated by sampling (once) :
∂ωk φ̃k∗ =  − 1√ Vk E ∼N (0,1)K ( −ωk∗+ √ Vk∗−ωk√ Vk ) k 6= k∗∑
k′ 6=k∗ 1√ Vk′
E ∼N (0,1)K ( −ωk∗+
√ Vk∗−ωk′√ Vk′
) k = k∗
(129)
where we have defined:
K(x) = N (x) H(x) =
√ 2/π
erfcx(x/2) (130)

heading: A.9.2 APPROACH 2: JENSEN AGAIN
text: A further simplification is obtained by applying Jensen inequality again to 128 but in the opposite direction, therefore we renounce to having a bound and look only for an approximation. We have the new effective free energy:
φ̃k∗ = ∑ k 6=k∗ logE ∼N (0,1)H ( −ωk ∗ + √ Vk∗ − ωk√ Vk ) (131)
= ∑ k 6=k∗ logH ( − ωk ∗ − ωk√ Vk + Vk∗ ) (132)
This gives, for k 6= k∗:
∂ωk φ̃k∗ =  − 1√ Vk+Vk∗ K ( − ωk∗−ωk√ Vk+Vk∗ ) k 6= k∗∑
k′ 6=k∗ 1√ Vk′+Vk∗ K ( − ωk∗−ωk′√ Vk′+Vk∗ ) k = k∗
(133)
Notice that ∂ωk∗ φ̃k∗ = − ∑ k 6=k∗ ∂ωk φ̃k∗ . In last formulas we used the definition of K in Eq. 130.
We show in Fig. 4 the negligible difference between the two ArgMax versions when using BP on the layers before the last one (which performs only the ArgMax).

heading: B EXPERIMENTAL DETAILS

heading: B.1 HYPER-PARAMETERS OF THE BP-BASED SCHEME
text: We include here a complete list of the hyper-parameters present in the BP-based algorithms. Notice that, like in the SGD type of algorithms, many of them can be fixed or it is possible to find a prescription for their value that works in most cases. However, we expect future research to find even more effective values of the hyper-parameters, in the same way it has been done for SGD. These hyper-parameters are: the mini-batch size bs; the parameter ρ (that has to be tuned similarly to the learning rate in SGD); the damping parameter α (that performs a running smoothing on the BP fields along the dynamics by adding a fraction of the field at the previous iteration, see Eqs. (134, 135)); the initialization coefficient that we use to to sample the parameters of our prior distribution qθ(W) according to θ`,t=0ki ∼ N (0, 1). Different choices of correspond to different initial distribution of the weights’ magnetization m`ki = tanh(θ ` ki), as is shown in Fig. 5); the number of internal steps of reinforcement τmax and the associated intensity of the internal reinforcement r. The performances of the BP-based algorithms are robust in a reasonable range of these hyper-parameters. A more principled choice of a good initialization condition could be made by adapting the technique from Stamatescu et al. (2020).
Notice that among these parameters, the BP dynamics at each layer is mostly sensitive to ρ and α, so that in general we consider them layer-dependent. See Sec. B.8 for details on the effect of these parameters on the learning dynamics and on layer polarization (i.e. how the BP dynamics tends to bias the weights towards a single point-wise configuration with high probability). Unless otherwise stated we fix some of the hyper-parameters, in particular: bs = 128 (results are consistent with other values of the batch-size, from bs = 1 up to bs = 1024 in our experiments), = 1.0, τmax = 1, r = 0.

heading: B.2 DAMPING SCHEME FOR THE MESSAGE PASSING
text: We use a damping parameter α ∈ (0, 1) to stabilize the training, changing the updated rule for the weights’ means as follows
m̃`,τki =∂Hψ(H `,τ−1 ki , G `,τ−1 ki , θ ` ki) (134) m`,τki =αm `,τ−1 ki + (1− α) m̃ `,τ ki (135)

heading: B.3 ARCHITECTURES
text: In the experiments in which we vary the architecture (see Sec. 4.1), all simulations of the BP-based algorithms use a number of internal reinforcement iterations τmax = 1. Learning is performed on the totality of the training dataset, the batch-size is bs = 128, the initialization coefficient is = 1.0.
For all architectures and all BP approximations, we use α = 0.8 for each layer, apart for the 501-501-501 MLP in which we use α = (0.1, 0.1, 0.1, 0.9). Concerning the parameter ρ, we use ρ = 0.9 on the last layer for all architectures and BP approximations. On the other layers we use: for the 101-101 and the 501-501 MLPs, ρ = 1.0001 for all BP approximations; for the 101- 101-101 MLP, ρ = 1.0 for BP and AMP while ρ = 1.001 for MF; for the 501-501-501 MLP ρ = 1.0001 for all BP approximations. For the BinaryNet simulations, the learning rate is lr = 10.0 for all MLP architectures, giving the better performance among the learning rates we have tested, lr = 100, 10, 1, 0.1, 0.001.
We notice that while we need some tuning of the hyper-parameters to reach the performances of BinaryNet, it is possible to fix them across datasets and architectures (e.g. ρ = 1 and α = 0.8 on each layer) without in general losing more than 20% (relative) of the generalization performances, demonstrating that the BP-based algorithms are effective for learning also with minimal hyperparameter tuning.
The experiments on the Bayesian error are performed on a MLP with 2 hidden layers of 101 units on the MNIST dataset (binary classification). Learning is performed on the totality of the training dataset, the batch-size is bs = 128, the initialization coefficient is = 1.0. In order to find the pointwise configurations we use α = 0.8 on each layer and ρ = (1.0001, 1.0001, 0.9), while to find the Bayesian ones we use α = 0.8 on each layer and ρ = (0.9999, 0.9999, 0.9) (these value prevent an excessive polarization of the network towards a particular pointwise configurations).
For the continual learning task (see Sec. 4.5) we fixed ρ = 1 and α = 0.8 on each layer as we empirically observed that polarizing the last layer helps mitigating the forgetting while leaving the single-task performances almost unchanged.
In Fig. 6 we report training curves on architectures different from the ones reported in the main paper.
B.4 VARYING THE DATASET
When varying the dataset (see Sec. 4.3), all simulation of the BP-based algorithms use a number of internal reinforcement iterations τmax = 1. Learning is performed on the totality of the training dataset, the batch-size is bs = 128, the initialization coefficient is = 1.0. For all datasets (MNIST (2 classes), FashionMNIST (2 classes), CIFAR-10 (2 classes), MNIST, FashionMNIST, CIFAR-10) and all algorithms (BP, AMP, MF) we use ρ = (1.0001, 1.0001, 0.9) and α = 0.8 for each layer. Using in the first layers values of ρ = 1 + with ≥ 0 and sufficiently small typically leads to good results.
For the BinaryNet simulations, the learning rate is lr = 10.0 (both for binary classification and multi-class classification), giving the better performance among the learning rates we have tested, lr = 100, 10, 1, 0.1, 0.001. In Tab. 2 we report the final train errors obtained on the different datasets.

heading: B.5 LOCAL ENERGY
text: We adapt the notion of flatness used in (Jiang et al., 2020; Pittorino et al., 2021), that we call local energy, to configurations with binary weights. Given a weight configuration w ∈ {±1}N , we define the local energy δEtrain(w, p) as the average difference in training error Etrain(w) when perturbingw by flipping a random fraction p of its elements:
δEtrain(w, p) = Ez Etrain(w z)− Etrain(w), (136)
where denotes the Hadamard (element-wise) product and the expectation is over i.i.d. entries for z equal to −1 with probability p and to +1 with probability 1− p. We report the resulting local energy profiles (in a range [0, pmax]) in Fig. 7 right panel for BP and BinaryNet. The relative error grows
slowly when perturbing the trained configurations (notice the convexity of the curves). This shows that both BP-based and SGD-based algorithms find configurations that lie in relatively flat minima in the energy landscape. The same qualitative phenomenon holds for different architectures and datasets.

heading: B.6 SGD IMPLEMENTATION (BINARYNET)
text: We compare the BP-based algorithms with SGD training for neural networks with binary weights and activations as introduced in BinaryNet (Hubara et al., 2016). This procedure consists in keeping a continuous version of the parameters w which is updated with the SGD rule, with the gradient calculated on the binarized configuration wb = sign(w). At inference time the forward pass is calculated with the parameters wb. The backward pass with binary activations is performed with the so called straight-through estimator.
Our implementation presents some differences with respect to the original proposal of the algorithm in (Hubara et al., 2016), in order to keep the comparison as fair as possible with the BP-based algorithms, in particular for what concerns the number of parameters. We do not use biases nor batch normalization layers, therefore in order to keep the pre-activations of each hidden layer normalized we rescale them by 1√
N where N is the size of the previous layer (or the input size in the case of the
pre-activations afferent to the first hidden layer). The standard SGD update rule is applied (instead of Adam), and we use the binary cross-entropy loss. Clipping of the continuous configuration w in [−1, 1] is applied. We use Xavier initialization (Glorot & Bengio, 2010) for the continuous weights. In Fig.2 of the main paper, we apply the Adam optimization rule, noticing that it performs slightly better in train and test generalization performance compared to the pure SGD one.

heading: B.7 EBP IMPLEMENTATION
text: Expectation Back Propagation (EBP) Soudry et al. (2014b) is parameter-free Bayesian algorithm that uses a mean-field (MF) approximation (fully factorized form for the posterior) in an online environment to estimate the Bayesian posterior distribution after the arrival of a new data point. The main differences between EBP and our approach relies in the approximation for the posterior distribution. Moreover we explicitly base the estimation of the marginals on the local high entropy structure. The fact that EBP works has no clear explanation: certainly it cannot be that the MF assumption holds for multi-layer neural networks. Still, it’s certainly very interesting that it works. We argue that it might work precisely by virtue of the existence of high local entropy minima and
expect it to give similar performance to the MF case of our algorithm. The online iteration could in fact be seen as way of implementing a reinforcement.
We implemented the EBP code along the lines of the original matlab implementation (https://github.com/ExpectationBackpropagation/EBP_Matlab_Code). In order to perform a fair comparison we removed the biases both in the binary and continuous weights versions. It is worth noticing that we faced numerical issues in training with a moderate to big batchsize All the experiments were consequently limited to a batchsize of 10 patterns

heading: B.8 UNIT POLARIZATION AND OVERLAPS
text: We define the self-overlap or polarization of a given unit a as qa0 = 1 N ∑ i(w a i )
2, where N is the number of parameters of the unit and {wai }Ni=1 its weights. It quantifies how much the unit is polarized towards a unique point-wise binary configuration (qa0 = 1 corresponding to high confidence in a given configurations while qa0 = 0 to low). The overlap between two units a and b (considered in the same layer) is qab = 1N ∑ wai w b i . The number of parameters N is the same for units belonging to the
same fully connected layer. We denote by 〈q0〉 = 1Nout ∑Nout a=1 q a 0 and 〈qab〉 = 1Nout ∑Nout a<b qab the mean polarization and mean overlap in a given layer (where Nout is the number of units in the layer).
The parameters ρ and α govern the dynamical evolution of the polarization of each layer during training. A value ρ ' 1 has the effect to progressively increase the units polarization during training, while ρ < 1 disfavours it. The damping α which takes values in [0, 1] has the effect to slow the dynamics by a smoothing process (the intensity of which depends on the value of α), generically favoring convergence. Given the nature of the updates in Algorithm 1, each layer presents its own dynamics given the values of ρ` and α` at layer `, that in general can differ from each other.
We find that it is is beneficial to control the polarization layer-per-layer, see Fig. 8 for the corresponding typical behavior of the mean polarization and the mean overlaps during training. Empirically, we have found that (as we could expect) when training is successful the layers polarize progressively towards q0 = 1, i.e. towards a precise point-wise solution, while the overlaps between units in each hidden layer are such that qab 1 (indicating low symmetry between intra-layer units, as expected for a non-redundant solution). To this aim, in most cases α` can be the same for each layer, while tuning ρ` for each layer allows to find better generalization performances in some cases (but is not strictly necessary for learning).
In particular, it is possible to use the same value ρ` for each layer before the last one (` < L where L is the number of layers in the network), while we have found that the last layer tends to
polarize immediately during the dynamics (probably due to its proximity to the output constraints). Empirically, it is usually beneficial for learning that this layer does not or only slightly polarize, i.e. 〈q0〉 1 (this can be achieved by imposing ρL < 1). Learning is anyway possible even when the last layer polarizes towards 〈q0〉 = 1 along the dynamics, i.e. by choosing ρL sufficiently large. As a simple general prescription in most experiments we can fix α = 0.8 and ρL = 0.9, therefore leaving ρ`<L as the only hyper-parameter to be tuned, akin to the learning rate in SGD. Its value has to be very close to 1.0 (a value smaller than 1.0 tends to depolarize the layers, without focusing on a particular point-wise binary configuration, while a value greater than 1.0 tends to lead to numerical instabilities and parameters’ divergence).

heading: B.9 COMPUTATIONAL PERFORMANCE: VARYING BATCH-SIZE
text: In order to compare the time performances of the BP-based algorithms with our implementation of BinaryNet, we report in Fig. 9 the time in seconds taken by a single epoch of each algorithm in function of the batch-size, on a MLP of 2 layers of 501 units on Fashion-MNIST. We test both algorithms on a NVIDIA GeForce RTX 2080 Ti GPU. Multi-class and binary classification present a very similar time scaling with the batch-size, in both cases comparable with BinaryNet. Let us also notice that BP-based algorithms are able to reach generalization performances comparable to BinaryNet for all the values of the batch-size reported in this section.

Review summaries:
This paper introduces a belief-propagation message-passing training algorithm for multi-layer neural networks. This algorithm is adapted to mini-batch training and biases distributions toward high entropy solutions. Empirical results show that neural networks with discrete weights and activations trained with this algorithm achieve comparable performance the same networks trained with SGD (BinaryNet), and can make approximate Bayesian predictions that have higher accuracy than pointwise solutions.

This manuscript provides an interesting try on alterative training algorithms for deep neural networks, based on (approximate) message-passing algorithms based on the well-known belief propagation (BP) algorithm.  In particular, the binary neural network is considered and four algorithms (BP, three variants of BP, i.e., BPI, MF, AMP) are proposed within a unified Posterior-As-Prior Update framework. Experiments are conducted on standard supervised classification tasks and continual learning settings, which shows comparable performances as standard SGD based methods.   ========================== After rebuttal: =========================  I have read the authors' feedback (many thanks for the detailed point-to-point feedback) and other reviewers' comments and modified the score accordingly. Overall, the proposed scheme is interesting, though strictly speaking the results are not very advantageous (at least from its current results) compared to traditional ones, and some of the comparisons seem not very reasonable/fair.