{'fraction_non_alphanumeric': 0.04158293015754869, 'fraction_numerical': 0.013000868378613075, 'mean_word_length': 4.305515335000658, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nonlocal optimization of binary neural networks

We explore training Binary Neural Networks (BNNs) as a discrete variable inference problem over a factor graph. We study the behaviour of this conversion in an under-parameterized BNN setting and propose stochastic versions of Belief Propagation (BP) and Survey Propagation (SP) message passing algorithms to overcome the intractability of their current formulation. Compared to traditional gradient methods for BNNs, our results indicate that both stochastic BP and SP find better configurations of the parameters in the BNN.


Introduction
Belief propagation (BP) is a message passing algorithm used to perform efficient inference and marginalization on probabilistic graphical models over a large set of interacting variables that have a tree-structure graph. BP works efficiently by utilizing the distributive law (Koller et al., 2009) to reduce the exponential complexity of inference or marginalization. Loopy BP (LBP) is an iterative message updating procedure that can be applied to loopy graphs. LBP can converge to one of the many possible fixed points in a loopy graph (Yedidia et al., 2005). These fixed points are often satisfactory approximations to the true values of interest (Murphy et al., 2013). Survey propagation is another message passing algorithm extending LBP that aims to capture all the fixed points of the BP algorithm weighed by their corresponding Bethe partition function (Mezard and Montanari, 2009).
The message passing algorithms can be cast conveniently on factor graphs (Kschischang et al., 2001), which are bipartite graphs comprising of variables of interest and factors that capture the local relations among the variables. Due to intractability of high-degree factors, message passing algorithms have only been generally studied in cases where each factor is connected to a small subset of variables, with exceptions involving certain classes of factors (Srinivasa et al., 2017;Braunstein and Zecchina, 2006;Tarlow et al., 2010). For instance, SP has mostly been applied in constraint satisfaction problems (CSPs) (Braunstein et al., 2005) and recently in Ising problems (Srinivasa et al., 2016). Extending BP and SP to high degree factors can open up avenues to study some interesting problems such as BNNs from a different perspective. This provides several advantages with respect to gradient-based approaches. First, one can use an inherently nonlocal optimization method such as SP to obtain solutions that are not accessible to local methods. Moreover, these methods naturally accommodate inherently discrete loss functions. These approaches lend themselves naturally to a Bayesian setting, where distributions over parameters, rather than point estimates, are obtained. BNNs are of much interest since they address some of the vast memory and computing power demands of deep neural networks. They are suited to real-time applications where neural networks are embedded in an energy-constrained environment (Peters and Welling, 2018). Recent works on BNNs are gradient-based and use stochastic gradient descent (SGD) for training. Binarizing the weights can happen either after training (Louizos et al., 2017), or during training Peters and Welling, 2018;Achterhold et al., 2018).
The contributions of this work are as follows: In Sec. 3.1, we provide a probabilistic interpretation of BP and propose a stochastic extension that allows (arbitrary) factors with any degree of connectivity. In Sec. 3.2, we propose two extensions of SP with stochastic messages that extend SP to high degree factors and variables. We demonstrate in Sec. A that our proposed models are consistent with exact BP and SP in small BNNs. In Sec. 4.1, we study the performance of our models in under-and over-parameterized conditions and demonstrate that they are able to obtain better solutions than gradient descent (GD). Finally in Sec. 4.2, we examine the performance of our models in high degree factors including multi layered perceptrons (MLPs) and convolutional networks on MNIST dataset and show that they perform well compared with SGD under these conditions.

Background
Here, we define the general problem setting, review BP and SP and set the notation that is used throughout the following sections. Consider a function F (w) over w = {w 1 , . . . , w N } with |w i | = K, i.e., each w i has cardinality D (assuming one of D possible values). We adopt the convenient notation used in Srinivasa et al. (2016), where underlined variables represent sets. Marginalization of F (w) over a set of variables indexed by α ⊂ {1, . . . , N }, w α , involves a sum with D N −|α| terms (|α| is the cardinality of w α ) and is intractable.

Belief Propagation
When F (w) factorizes as a product of local "factors" f I (w I ), the distributive law (Bishop, 2006), ab + ac = a(b + c) can be used to solve the marginalization problem. The problem can be represented as messages passed on a factor graph (Kschischang et al., 2001) which is a bipartite graph comprising of factors, f I , and variables, w i , indexed by upper and lower letters, respectively. The computational gain is achieved by first applying a summation to the factors (instead of the whole F (w)) and then accumulating the results as products on variables. This "message-passing" algorithm leads to two types of messages, the variable to factor and factor to variable messages. The message, m I→i (w i ), passed from factor I to variable i is defined as: where the symbol ∂[.] denotes the set of nodes (variables or factors) adjacent to node [.], ∂I \ i → I represents all variables adjacent to factor I except i. The message from variable i to factor I is defined as: In order to make the equations less cluttered, the functional dependence of the messages on m ∂i\I→i and m ∂i\i→I are suppressed when no confusion arises. The belief of variable i is defined as which is equal to the true marginal in tree-structured factor graphs. BP works by choosing a variable as the root node in the graph, and subsequently passing messages between the variables and factors using Eqs. 1-2. When the factor graph is a tree, the computed marginal is exact. LBP can be used to perform approximate inference on loopy graphs by initializing the messages and iteratively applying Eqs. 1-2.

Survey Propagation
It is demonstrated by Heskes (2002) that when LBP converges (a stable fixed point) on a loopy factor graph, a local minimum of the Bethe approximation to the partition function is obtained. A factor graph with loops can have many BP fixed points. SP is a powerful extension of BP designed to keep track of an exponential number of fixed points in LBP (Ravanbakhsh and Greiner, 2014). This is achieved by assuming that the joint distribution represented by the factor graph can be approximated by a sum over the joint distributions obtained by fixed points of LBP. Therefore, marginalization in SP involves a sum over all corresponding variables for a given fixed point, and eventually, summing over all the fixed points. The conditions for validity of SP approximation correspond to a phase in replica symmetry breaking (Mezard and Montanari, 2009) known as the clustering phase which is not easy to check. Srinivasa et al. (2016) demonstrated that regardless of the validity of assumptions, SP can operate on a wide range of problems and obtain solutions that are not accessible to BP. Similar to BP, SP comprises of two sets of messages. The variable to factor message is defined as (Srinivasa et al., 2016): Here, the integration is over the space of all possible BP factor (∂i \ I) to variable (i) message functions m ∂i\I→i . M i→I (m i→I = m 0 ) can be interpreted as the message that variable i sends to factor I signaling that the corresponding variable to factor BP message is equal to m 0 . The indicator function I[.] constrains the integration region to domains that respect the corresponding BP update. Note that m i→I (.; m ∂i\I→i )(∅) wi m i→I (w i ; m ∂i\I→i ). Similarly, the factor to variable message is defined as (with similar notation and definitions): In SP, these two types of messages are updated iteratively until convergence and eventually used to calculate the marginals (See Sec. 3.2).
3 Proposed approach

Stochastic BP (SBP)
The factor to variable message update in BP, m I→i (w i ), (i.e., Eq. 1) involves a sum over w \i . However, one only needs to consider the variables w I , since f I (w I ) does not depend on the other variables w \i\I . Nonetheless, this sum involves D |w I | terms, which is only tractable for low-degree factors. In this section, we propose a stochastic formulation of BP that allows high degree factors.
Notice that the variable to factor messages that appear in Eq. 1 are all positive numbers, therefore, we can convert j∈∂I\i m j→I (w j ) to a joint probability distribution p(w ∂I\i ) j∈∂I\i p j→I (w j ) where p j→I (w j ) m j→I (wj ) m j→I (∅) and m j→I (∅) wj m j→I (w j ). This transformation is justified since we are only interested in the ratios of m I→i (w i ) for the possible D values of w i and this transformation multiplies all the possible values of the message by the same quantity. Therefore, Eq. 1 can be rewritten as: wherep Our proposed approach is to regard Eq. 5 as expectation value of f I (w I ) with respect to the joint probability distribution p(w ∂I\i ). The value of the message can be estimated using a Monte Carlo (MC) approach, i.e., drawing samples from the joint distribution: where m indexes over the L draws from the joint distribution. We call this model stochastic BP (SBP). It is not related to (Noorshams and Wainwright, 2013), where a stochastic approach is employed to deal with high values of D. Algorithm 1 illustrates this approach when the variables are binary (w i ∈ {0, 1}). This corresponds to a joint distribution over a factorial Bernoulli distribution. The algorithm loops over all variables i adjacent to a factor I (lines 1-11). It also loops over L MC samples drawn from the joint probability distribution defined by the variable to factor messages (lines 2-7). The ith elements, w i of the drawn samples are set to 1 and 0 at lines 4 and 6, respectively. These values are used to obtain MC estimates of the unnormalized probabilities of sending a 1 or a 0 from I to i in lines 8-9. The normalized probability of sending a 1 is calculated in line 10.
Algorithm 1: Stochastic BP factor to variable update subroutine for binary variables:

end
3.2 Stochastic SP (S 3 P and S 4 P) The integrals in the SP message updates are converted into sums by keeping histograms of BP messages. With K bins, the sums involve K |∂I|−1 and K |∂i|−1 terms, for variable to factor and factor to variable messages, respectively. These sums are both intractable for high-degree nodes. The exponential complexity of the variable to factor messages can be reduced to O(K 2 |∂i|) by performing FFT, multiplying the signals, and then doing an inverse FFT, see (Srinivasa et al., 2016) for details. This process involves padding the original messages from size K to ≈ |∂i|K in order to avoid aliasing. Since |∂i| represents the number of examples in a dataset (see Sec. 4), it is not practical for large datasets. Here, we provide probabilistic interpretations of both these messages to obtain stochastic estimates of the messages in linear time.
Consider the SP factor to variable update, Eq. 4. By definingS j→I (p j→I ) ∝ M j→I (p j→I ) with S j→I (p j→I ) S j→I (p j→I ) p j→IS j→I (p j→I ) , we can regard the factor to variable update as an expectation under the joint distribution S(p ∂I\i ) j∈∂I\i S j→I (p j→I ). Note that each S j→I (p j→I ) is a categorical distribution (with cardinality K). The resulting factor to variable message is: where p k ∈ [0, 1] is one of the possible K values that p I→i can take. p I→i (.; p ∂I\i→I )(∅) is estimated stochastically using Algorithm 1. Algorithm 2 shows how this approach works. The algorithm receives variable to factor surveys (SP messages) as input {S i→I (p i→I )} {i∈∂I} . The MC estimate is done over L sp samples (lines 1-9). A variable to factor probability distribution is sampled from the surveys incident on factor I. This probability distribution is fed into the StochasticBPF2Var subroutine to obtain the normalized and unnomarlized probabilities (line 3). These unnormalized probabilities are used to calculate the local partition function for each variable (line 5), whereas the normalized probability is used to find the bin k that corresponds to the given probability (line 6), where bin(p k ) k. The message corresponding to this bin is subsequently updated (line 7). Finally, the normalized messages are calculated (lines 10-14). A similar approach can be used to update the variable to factor messages (See Appendix B). We call our model S 4 P (the first 3 S's correspond to the three sources of stochasticity from the factor to variable, variable to factor and BP messages, while the last S stands for survey) as opposed to S 3 P that uses the FFT trick to calculate the variable to factor messages.
Since SP is a generalization of the BP algorithm, we call a combination of SBP and S 4 P the Stochastic Nonlocal Message Passing (SNMP) model. SNMP starts by running the faster BP first, and if all the constraints are not satisfied, S 4 P is called.
Algorithm 2: SP Stochastic factor to variable subroutine.

Problem set-up
Consider a training set of binary vectors x and labels, y, comprising of M examples. We assume binary labels, y I ∈ {0, 1}, for simplicity in this work, though any number of labels can also be studied using this approach. Assume a binary neural network with weights w that predicts label g w (x I ) given input x I . g w (.) can be any function including MLPs with any type of activation. We assume that the last layer of MLP is a sign function. The training of the MLP can be posed as a CSP problem using a (high degree) factor graph containing M factors and N |w| variables: F (w) = M I=1 f I (w) and f I (w) I g w (x I ) = y I . We have considered the general case where each factor f I depends on all the variables, i.e., w I = w. With this set-up, BP and SP algorithms can be employed to calculate the marginals of the weights of the neural networks that satisfy all the instances on the training set. We relax the strict requirement of satisfying all the training examples by adopting a Max-Sat strategy (Chieu and Lee, 2009): in sum-product forms of BP and SP, F (w) is zero even when one of the factors is not satisfied. This makes the algorithm attempt to satisfy all the instances and not directly discriminate between cases that cannot satisfy one or more of the constraints. We relax this by writing Here, β is the inverse temperature for a Boltzmann distributions in the language of statistical mechanics (assuming Boltzmann constant is unity). A high value of β penalizes the unsatisfied factors more severely.

Experiments
We first demonstrate in Fig. 1 that SBP and S 4 P yield comparable results to the BP and S 3 P, respectively on a whole range of different conditions (see Appendix. A for details). We do not compare with exact SP since the factor to variable messages are intractable even for small problems described below. Our experiments are implemented in pytorch (Paszke et al., 2017). We use broadcasting operations on the GPU to parallelize all operations. Performance of S 4 P against S 3 P: for larger values of samples S 4 P is able to reach higher accuracies, due to accumulative binning errors that could arise in S 3 P (see Appendix. A for details).

Nonlocal solver under different regimes
Here, we investigate the utility of message passing algorithms for classification tasks under carefully designed conditions. These conditions are chosen so that the difference between local and nonlocal solvers can be highlighted. Consider what we call in our experiments the glass dataset where x and y comprise of independent fair Bernoulli units. We are interested in neural nets under over-parametrized and under-parametrized conditions, and the region between them. The parameter α M N , i.e., the ratio of the number of training examples to the number of variables is a measure of how over-or under-parameterized the system is (Mezard and Montanari, 2009) for this dataset. When the system is over-parameterized (small values of α), many solutions exist. For an under-parameterized system at large values of α, it is impossible to satisfy all the constraints. The phase transition between the two is of great interest, since the solutions are difficult to find in this region; a solver that undergoes phase transition at a higher value of α is therefore more powerful (Chieu and Lee, 2009). It is known that neural networks with full precision weights with some mild conditions on the number of hidden units and activations are easily able to learn any random assignment of labels to inputs (Zhang et al., 2017). Therefore, in order to clearly see the over-and under-parameterization, we use a single one-layer MLP with a sign function output. Fig. 4(left) depicts the performance of different algorithms on the glass dataset. Note that SNMP's accuracy drops at higher values of α compared with the local methods SBP and GD. The numbers on the vertical axis are the average accuracies across 20 lines. See Appendix C.2 for details.
As another benchmark, and in contrast to the glass dataset with no structure in its elements, we study the digits dataset (Alpaydin and Alimoglu, 1998) comprising 250 images with dimensionality of 64. The results are portrayed in Fig. 4(right). Note that due to simplicity and structure of this dataset, a linear model is not under-parameterized and the difference between the global and local methods are less pronounced.
As an intermediary between these two datasets, we designed the stained glass dataset, where the inputs are sampled from fair Bernoulli distributions and the output is generated by passing the input through an MLP (with 1 or 2 hidden layers and sign non-linearity) with random binary weights. The weights are random binary variables, and the labels, y, have a functional relationship with x; a relationship that is not necessarily captured by a single-layered MLP. Fig. 4(middle) compares the performance of different models under this dataset. Interestingly, the gap between the local and nonlocal methods is in between the two extreme cases, suggesting that nonlocal binary solvers can be of more utility in underparameterized settings.
Amongst possible approaches to using (S)GD for optimization of BNNs, we have adopted the method of Courbariaux and Bengio (2016). This approach does not rely on the central limit theorem (CLT) as do methods by Peters and Welling (2018); Shayer et al. (2018). CLT breaks when the number of hidden units is not large enough which is the case in a one-layer network. Moreover, Courbariaux and Bengio (2016) use the full precision of weights in back propagation, and only use binary weights during the forward pass resulting in competitive performance.

Benefiting from non-locality
One of the specific advantages of using SP over a local solver like BP or GD is its ability to efficiently consider a large number of candidate solutions. Because SP is able to track and aggregate information from potentially an exponential number of BP fixed points (Ravanbakhsh and Greiner, 2014), it can often find solutions that are inaccessible to a single SBP or GD run. We demonstrate this in Figure  3, which shows the evolution of a given S 4 P marginal over iterations of the algorithm. Each box represents a model, initialized in the same way and trained on the same dataset, but by varying the number of samples used in stochastic SP, we can control the ability of the algorithm to explore the solution space. The extreme case, L sp = 1, is also presented, and contrasted with SBP. In the case of BP, we see the beliefs hop around until they eventually converge on a single solution which achieves sub-optimal accuracy. As we increase the sample size, more and more mixing is possible, which translates to more accurate configurations that were previously inaccessible to the model.  BP is presented for contrast. Each box represents the time evolution of the SP marginal for a given variable. While BP can only track (and converge to) one candidate solution, SP is free to explore and propagate multiple candidates non-locally. As L sp increases S 4 P can find better configurations which are inaccessible to BP. Best viewed in colour.

MNIST dataset
We study the performance of our SNMP against gradient-based methods using a binary classification task on MNIST dataset. We have chosen linear, MLP and convolutional networks. See Appendix C.3 for details. The phase transition diagrams for both the test and training sets are plotted for the linear, MLP and convolutional networks in Fig. 4. In order to scale our vectorized implementation to a larger number of MC samples, we use a mini-batch version of message passing as delineated in Appendix. C.3. This allows running the algorithm on all examples of MNIST (from the given classes).  Figure 5: Training curves on full MNIST dataset for multi-layer convolutional network trained with Max-Sat, using a mini-batch approach. We observe that higher β leads to a higher accuracy score, though noticeably β ∞ is not the optimal value. Shown for comparison is SGD. The figure on the left is the training accuracy and on the right is the test accuracy.

Conclusions and future work
We extended the framework of message passing to arbitrary factors with high degrees and used this approach to obtain nonlocal solvers for binary neural networks. Our studies reveal that using powerful nonlocal solvers can be beneficial in underparameterized systems. These systems are interesting in energy-constrained environments where vast computational resources are not available. We showed that our models perform competitively with gradient based methods and are able to access parts of the solution space that is not directly accessible to local methods. Application and extension of these models to large scale datasets is a prime direction for our future studies.

Supplementary Material for Nonlocal optimization of binary neural networks A Consistency between stochastic and exact models
We gauge the accuracy of SBP and S 4 P against BP and S 3 P on a large number of problems and range of parameters for which BP and S 3 P are tractable. More specifically, we fix M = 2 N , and sample all combinations from a grid formed by: • N ∈ {4, 6, 8, 10} • γ ∈ {0.2, 0.5, 0.8} • L BP ∈ {5, 8, 10, 12, 15, 20, 25, 50} SBP and BP are run on the resulting problems. A comparison between the accuracies can be seen in Fig. 1(left), where we find nearly identical agreement even when a small number of samples are used. It can also be seen that SBP is relatively insensitive to the choice of L BP , i.e., the number of MC samples.

B Stochastic factor to variable update in SP
Similar to the derivations in the main text, the variable to factor message can be written as: where p k ∈ [0, 1] is one of the possible K values that p i→I can take. p I→i (.; p ∂I\i→I )(∅) is calculated directly using Eq. 2. Algorithm 3 shows how this approach works. The algorithm receives variable to Factor surveys (SP messages) as input {S i→I (p i→I )} {i∈∂I} . The MC estimate is done over L sp samples (lines 1-9). A variable to factor probability distribution is sampled from the surveys incident on variable i. This probability distribution is fed into the BPVar2F subroutine (Eq. 2) to obtain the normalized and unnomarlized probabilities (line 3). These unnormalized probabilities are used to calculate the local partition function for each variable (line 5), whereas the normalized probability is used to find the bin k that corresponds to the given probability (line 6), where bin(p k ) k. The message corresponding to this bin is subsequently updated (line 7). Finally, the normalized messages are calculated (lines 10-14).

C Experimental and implementation details C.1 General considerations when tuning hyperparameters
Here, we summarize the different hyperparameters used in our solver, and their recommended settings.  • L SP : the number of samples used for S 3 P and S 4 P. As shown in Fig. 3, this parameter is crucial for proper exploration of the space, and using higher L SP will typically improve the final accuracy. From our own tuning, we empirically find that one wants at least L SP > 10L BP . If the memory space is available, L SP is one of the first things to tune. In all our experiments we have at least L SP = 100, though we go up to 500 when we can afford it (as done for Fig. 3).
• K: the number of bins used for SP. A low value results in insufficient resolution to track the different BP messages. However, adjusting K offers little after a certain threshold.
• γ: the damping factor. We use damping in our updates to make the system more stable (Murphy et al., 2013). If the previous value of a message (variable to factor or otherwise) is given by µ t , then with damping factor γ after computing updates µ, the new message becomes: µ t+1 = (1 − γ)µ t + µ. When L SP is high, γ can afford to be higher, and converges more quickly. For SBP we typically set γ = 0.2, while for S 4 P we use γ = 0.8.
• β: the inverse temperature for Max-Sat. Higher values of β have a positive impact on the variance encountered during training, especially when combined with tuning γ. Fig. 5 provides a good illustration of the impact of β on training.

C.2 Details for experiments in Sec. 4.1
The message passing algorithms converge much faster (usually less than 5 epochs) than gradient based methods. However, we run all the algorithms for a higher number of epochs (20) so that SGD can reach its solution. We fine-tune the learning rate of SGD for optimum performance. The best performance was obtained with a learning rate of 0.1.

C.3 Details for experiments in Sec 4.2
We test the solver's ability to handle models beyond simple linear classifiers, and test on a larger dataset, namely MNIST. The linear model and the MLP operate on flattened images. In order to reduce the number of variables required for this, we use a max-pooling operation, reducing the images to size (14,14), and thus the input dimension to 14 × 14 = 196. The linear model thus consists of 196 parameters, while the MLP, with a single hidden layer of size 2, has 196 × 2 + 2 = 394 parameters.
Finally, the conv net consists of three layers: • A Convolutional layer with 3 (3x3) filters, using a stride of 2. The activation is a simple threshold y = (h > 0) • A second convolutional layer, with 2 (3x3) filters, also using a stride of 2, followed by threshold at 0.
• The outputs of the previous layer are flattened to form a vector of size (8), after which a linear layer with a sign function activation is used.
Taken together, the convolutional network has 89 parameters. Note that none of our models contain biases.
We use γ = 0.5 for the linear model and the MLP. For the convolutional network, we adapt γ according to the present accuracy of the model. At the beginning of training, when accuracy is low, γ should be high to explore the space quickly. As the model becomes more accurate, however, γ can be reduced to limit the variance. As such, every 20 iterations we compute the accuracy, and set γ = (1 − accuracy). As Sec.C.4 describes, our implementation vectorizes the message computation, and performs all MC estimations in one pass. The tradeoff is that this requires us to multiply the size of those messages in memory by L SP . Thus, rather than using a smaller sample size for our MC estimation, we draw inspiration from minibatch optimization and consider ways to update the messages using a smaller subset of the factors. At a given iteration, we select a subset of the factors, and only update the factor to variable messages from those factors. We then immediately update all of the variable to factor messages, instead of looping over the remaining factors. Fig. 5 shows the training history of a three layer convolutional neural network, using the Max-Sat strategy described in Sec. 3.3. We explore β in the range {0.1, 1, 10, 100, ∞}, where ∞ corresponds to regular SBP.

C.4 General considerations for vectorized implementations
For BP and SBP, the variable to factor and factor to variable messages can be represented as matrices with shape [M, N ], and similarly, for S 3 P and S 4 P the messages have shape [M, N, K]. The main computational demands of the algorithm are in the generation, and subsequent evaluation, of candidate configurations for the variables.
For generating configurations, the main use of memory is in sampling large multi-dimensional matrices. For instance, to compute the message updates in SBP described in Alg. 1, each factor requires sampling a matrix with shape [L BP , 2N, N ], which grows quadratically with the number of variables. We can generate the full set of samples in one pass by uniformly sampling a 4D matrix of size [M, L BP , 2N, N ], and then performing broadcasted element-wise comparison with the variable-to-factor messages along the second and third dimensions of the matrix.
For the factor-to-variable messages in Alg. 2, before we can sample the candidate configurations, we must first generate samples from the categorical distribution over the K bins described by the variable to factor messages. We may generate these by sampling uniform matrices with shape [L SP , M, N ], and performing broadcasted comparisons with the cumulative distribution over the K bins. Having obtained these sample bin indices, we can look up the corresponding probabilities, which we can use to generate the sample configurations. This proceeds in the same way as SBP, except we now sample 5D matrices of shape [M, L SP , L BP , 2N, N ].
The configurations, as sampled, consist of flattened and concatenated values for the parameters of the BNN. We load and reshape contiguous subsets of these parameters according to the architecture of the model. For simple linear layers, matrix mulplication can be performed by broadcasting element wise multiplication, followed by summation. The evaluation of convolutional layers can be performed efficiently using grouped convolution. If the convolutional layer requires filters of shape [C out , C in , K W , K H ], the corresponding subset of the configuration matrix can be reshaped to [M, L SP , L BP , 2N, C out , C in , K W , K H ]. To leverage the efficient grouped convolutions present in pytorch, we reshape these configurations into a matrix of shape The inputs to the layer are reshaped to [1, B × C in , W, H]. We then perform convolution using groups of size B (and subsequently unpack the reshaped outputs).
With the batch version, we seek a tradeoff between memory requirements and the amount of computation. We sample a subset of the factors of size b, and use the same procedure above to compute the messages to all variables. The matrices sampled in this way posses the same dimensionality.
We gauge the accuracy of SBP and S 4 P against BP and S 3 P on a large number of problems and range of parameters for which BP and S 3 P are tractable. More specifically, we fix M = 2 N , and sample all combinations from a grid formed by: , 8, 10, 12, 15, 20, 25, 50} SBP and BP are run on the resulting problems. A comparison between the accuracies can be seen in Fig. A.1(left), where we find nearly identical agreement even when a small number of samples are used. It can also be seen that SBP is relatively insensitive to the choice of L BP , i.e., the number of MC samples.
In a similar way, to compare S 3 P and S 4 P, we fix N = 20 and sample all combinations from a grid of: A comparison between these accuracies is also present on the right side of Fig. A.1. Two important points are noted: first, higher numbers of MC samples lead to better S 4 P performance. Second, quite interestingly, S 4 P usually outperforms S 3 P, due to the binning errors that accumulate if all the the exponential combinations of bins are considered.

B Stochastic factor to variable update in SP
Similar to the derivations in the main text, the variable to factor message can be written as: where p k ∈ [0, 1] is one of the possible K values that p i→I can take. p I→i (.; p ∂I\i→I )(∅) is calculated directly using Eq. ??. Algorithm 1 shows how this approach works. The algorithm receives variable to Factor surveys (SP messages) as input {S i→I (p i→I )} {i∈∂I} . The MC estimate is done over L sp samples (lines 1-9). A variable to factor probability distribution is sampled from the surveys incident on variable i. This probability distribution is fed into the BPVar2F subroutine (Eq. ??) to obtain the normalized and unnomarlized probabilities (line 3). These unnormalized probabilities are used to calculate the local partition function for each variable (line 5), whereas the normalized probability is used to find the bin k that corresponds to the given probability (line 6), where bin(p k ) k. The message corresponding to this bin is subsequently updated (line 7). Finally, the normalized messages are calculated (lines 10-14).  Here, we summarize the different hyperparameters used in our solver, and their recommended settings.
• L BP : the number of samples used for SBP. As shown on the left of Fig. A.1, one can get away with typically very small number of samples. In all our experiments, we set L BP = 5 • L SP : the number of samples used for S 3 P and S 4 P . As shown in Fig. ??, this parameter is crucial for proper exploration of the space, and using higher L SP will typically improve the final accuracy. From our own tuning, we empirically find that one wants at least L SP > 10L BP . If the memory space is available, L SP is one of the first things to tune. In all our experiments we have at least L SP = 100, though we go up to 500 when we can afford it (as done for Fig. ??).
• K: the number of bins used for SP. A low value results in insufficient resolution to track the different BP messages. However, adjusting K offers little after a certain threshold.
• γ: the dampening factor. Helps to reduce the variance, especially when adaptive methods are used, as in Sec. C.3. When L SP is high, γ can afford to be higher, and converges more quickly.
• β: the inverse temperature for Max-Sat. Higher values of β have a positive impact on the variance encountered during training, especially when combined with tuning γ. Fig. ?? provides a good illustration of the impact of β on training.

C.2 Details for experiments in Sec. ??
The message passing algorithms converge much faster (usually less than 5 epochs) than gradient based methods. However, we run all the algorithms for a higher number of epochs (20) so that SGD can reach its solution. We fine-tune the learning rate of SGD for optimum performance. The best performance was obtained with a learning rate of 0.1. For Fig. ??, we set N = 10, M = 20, and ran for 20 iterations with γ = 0.5. K was set to 31 for ease of visualization.

C.3 Details for experiments in Sec ??
As Sec.C.4 describes, our implementation vectorizes the message computation, and performs all MC estimations in one pass. The tradeoff is that this requires us to multiply the size of those messages in memory by L sp . Thus, rather than using a smaller sample size for our MC estimation, we draw inspiration from minibatch optimization and consider ways to update the messages using a smaller subset of the factors. At a given iteration, we select a subset of the factors, and update only the messages sent out by those factors. When then immediately update all of the variable to factor messages, instead of looping over the remaining factors. We use damping in our updates to make the system more stable (?). If the previous value of a message (variable to factor or otherwise) is given by µ t , then with damping factor γ after computing updates µ, the new message becomes: µ t+1 = (1 − γ)µ t + µ. In order to make the runs more memory efficient, we use a max-pooling operation, reducing the images to size (14,14). The linear model thus consists of 14 × 14 = 196 parameters, while the MLP, with a single hidden layer of size 2, has 196 × 2 + 2 = 394 parameters. Finally, the conv net consists of three layers: • A Convolutional layer with 3 (3x3) filters, using a stride of 2. This is followed by thresholding at 0. • A second convolutional layer, with 2 (3x3) filters, also using a stride of 2, followed by threshold at 0. • The outputs of the previous layer are flattened to form a vector of size (8), after which a linear layer with a sign function activation is used.
Note that none of our models contain biases. Taken together, the convolutional network has 89 parameters.
We demonstrate this flexibility on a much larger subset of the MNIST training set, and in addition, test the solver's ability to handle models beyond simple linear classifiers. Fig. ?? shows the training history of a three layer convolutional neural network, using the Max-Sat strategy described in Sec.
??. We explore β in the range {0.1, 1, 10, 100, ∞}, where ∞ corresponds to regular SBP. Via tuning we find the best value to be β = 100. We compare this, again, with a full-precision Adam optimizer, which converges within to 90% after 120 epochs.

C.4 General considerations for vectorized implementations
For BP and SBP, the variable to factor and factor to variable messages can be represented as matrices with shape [M, N ], and similarly, for S 3 P and S 4 P the messages have shape [M, N, K]. The main computational demands of the algorithm are in the generation, and subsequent evaluation, of candidate configurations for the variables.
For generating configurations, the main use of memory is in sampling large multi-dimensional matrices. For instance, to compute the message updates in SBP described in Alg. ??, each factor requires sampling a matrix with shape [L BP , 2N, N ], which grows quadratically with the number of variables. We can generate the full set of samples in one pass by uniformly sampling a 4D matrix of size [M, L BP , 2N, N ], and then performing broadcasted element-wise comparison with the variable-to-factor messages along the second and third dimensions of the matrix.
For the factor-to-variable messages in Alg. ??, before we can sample the candidate configurations, we must first generate samples from the categorical distribution over the K bins described by the variable to factor messages. We may generate these by sampling uniform matrices with shape [L SP , M, N ], and performing broadcasted comparisons with the cumulative distribution over the K bins. Having obtained these sample bin indices, we can look up the corresponding probabilities, which we can use to generate the sample configurations. This proceeds in the same way as SBP, except we now sample 5D matrices of shape [M, L SP , L BP , 2N, N ].
The configurations, as sampled, consist of flattened and concatenated values for the parameters of the BNN. We load and reshape contiguous subsets of these parameters according to the architecture of the model. For simple linear layers, matrix mulplication can be performed by broadcasting element wise multiplication, followed by summation. The evaluation of convolutional layers can be performed efficiently using grouped convolution. If the convolutional layer requires filters of shape [C out , C in , K W , K H ], the corresponding subset of the configuration matrix can be reshaped to [M, L SP , L BP , 2N, C out , C in , K W , K H ]. To leverage the efficient grouped convolutions present in pytorch, we reshape these configurations into a matrix of shape [B × C out , C in , K W , K H ], where B = M × L SP × L BP × 2N . The inputs to the layer are reshaped to [1, B × C in , W, H]. We then perform convolution using groups of size B (and subsequently unpack the reshaped outputs).
With the batch version, we seek a tradeoff between memory requirements and the amount of computation. We sample a subset of the factors of size b, and use the same procedure above to compute the messages to all variables. The matrices sampled in this way posses the same dimensionality.