# Nonlocal optimization of binary neural networks

Amir Khoshaman
Borealis AI
&Giuseppe Castiglione11footnotemark: 1
Borealis AI
&Christopher Srinivasa
Borealis AI
Authors contributed equally to this work.

###### Abstract

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.

##  1 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 (Ullrich et al., 2017; 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.

##  2 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(\underline{w})$ over $\underline{w}=\\{w_{1},\dots,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(\underline{w})$ over a set of variables indexed by
$\alpha\subset\\{1,\dots,N\\}$, $\underline{w}_{\alpha}$, involves a sum with
$D^{N-|\alpha|}$ terms ($|\alpha|$ is the cardinality of
$\underline{w}_{\alpha}$) and is intractable.

###  2.1 Belief Propagation

When $F(\underline{w})$ factorizes as a product of local “factors”
$f_{I}(\underline{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, $\underline{f}_{I}$,
and variables, $\underline{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(\underline{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\rightarrow i}(w_{i})$, passed from factor $I$
to variable $i$ is defined as:

${m}_{I\rightarrow i}(w_{i};m_{\partial I\setminus i\rightarrow
I})\triangleq\sum_{\underline{w}_{\setminus
i}}f_{I}(\underline{w}_{I})\prod_{j\in\partial{I}\setminus{i}}m_{j\rightarrow
I}(w_{j}),$ (1)

where the symbol $\partial[.]$ denotes the set of nodes (variables or factors)
adjacent to node $[.]$, $\partial I\setminus i\rightarrow I$ represents all
variables adjacent to factor $I$ except $i$. The message from variable $i$ to
factor $I$ is defined as:

${m}_{i\rightarrow I}(w_{i};m_{\partial i\setminus I\rightarrow
i})\triangleq\prod_{J\in\partial{i}\setminus{I}}m_{J\rightarrow i}(w_{j}).$
(2)

In order to make the equations less cluttered, the functional dependence of
the messages on $m_{\partial i\setminus I\rightarrow i}$ and $m_{\partial
i\setminus i\rightarrow I}$ are suppressed when no confusion arises. The
belief of variable $i$ is defined as ${m}_{i}(w_{i};m_{\partial i\rightarrow
i})\triangleq\prod_{J\in\partial{i}}m_{J\rightarrow i}(w_{j})$ 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.

###  2.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):

$\displaystyle{M}_{i\rightarrow I}({m}_{i\rightarrow I}=m_{0})$
$\displaystyle\triangleq\int_{{\underline{m}^{\prime}}_{\partial i\setminus
I\rightarrow i}}\Bigg{[}\mathbb{I}\left[m_{0}={m}_{i\rightarrow
I}(.;m^{\prime}_{\partial i\setminus I\rightarrow i})\right]\times$
$\displaystyle{m}_{i\rightarrow I}(.;m^{\prime}_{\partial i\setminus
I\rightarrow i})(\emptyset)\prod_{J\in\partial i\setminus I}{M}_{J\rightarrow
I}({m^{\prime}}_{J\rightarrow i})\Bigg{]}.$ (3)

Here, the integration is over the space of all possible BP factor ($\partial
i\setminus I$) to variable ($i$) message functions
$\underline{m}^{\prime}_{\partial i\setminus I\rightarrow i}$.
${M}_{i\rightarrow I}({m}_{i\rightarrow 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
$\mathbb{I}[.]$ constrains the integration region to domains that respect the
corresponding BP update. Note that ${m}_{i\rightarrow
I}(.;m^{\prime}_{\partial i\setminus I\rightarrow
i})(\emptyset)\triangleq\sum_{w_{i}}{m}_{i\rightarrow
I}(w_{i};m^{\prime}_{\partial i\setminus I\rightarrow i})$. Similarly, the
factor to variable message is defined as (with similar notation and
definitions):

$\displaystyle{M}_{I\rightarrow i}({m}_{I\rightarrow i}=m_{0})$
$\displaystyle\triangleq\int_{{\underline{m}^{\prime}}_{\partial I\setminus
i\rightarrow I}}\Bigg{[}\mathbb{I}\left[m_{0}={m}_{I\rightarrow
i}(.;m^{\prime}_{\partial I\setminus i\rightarrow I})\right]\times$
$\displaystyle{m}_{I\rightarrow i}(.;m^{\prime}_{\partial I\setminus
i\rightarrow I})(\emptyset)\prod_{j\in\partial I\setminus i}{M}_{j\rightarrow
I}({m^{\prime}}_{j\rightarrow I})\Bigg{]}.$ (4)

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

###  3.1 Stochastic BP (SBP)

The factor to variable message update in BP, ${m}_{I\rightarrow i}(w_{i})$,
(_i.e._ , Eq. 1) involves a sum over $\underline{w}_{\setminus i}$. However,
one only needs to consider the variables $\underline{w}_{I}$, since
$f_{I}(\underline{w}_{I})$ does not depend on the other variables
$\underline{w}_{\setminus i\setminus I}$. Nonetheless, this sum involves
$D^{|\underline{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
$\prod_{j\in\partial{I}\setminus{i}}m_{j\rightarrow I}(w_{j})$ to a joint
probability distribution
$p({w}_{\partial{I}\setminus{i}})\triangleq\prod_{j\in\partial{I}\setminus{i}}{p}_{j\rightarrow
I}(w_{j})$ where ${p}_{j\rightarrow I}(w_{j})\triangleq\frac{{m}_{j\rightarrow
I}(w_{j})}{{m}_{j\rightarrow I}(\emptyset)}$ and ${m}_{j\rightarrow
I}(\emptyset)\triangleq\sum_{w_{j}}{m}_{j\rightarrow I}(w_{j})$. This
transformation is justified since we are only interested in the ratios of
${m}_{I\rightarrow 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:

$\tilde{p}_{I\rightarrow i}(w_{i};p_{\partial I\setminus i\rightarrow
I})=\sum_{\underline{w}_{\setminus
i}}f_{I}(\underline{w}_{I})\prod_{j\in\partial{I}\setminus{i}}{p}_{j\rightarrow
I}(w_{j}),$ (5)

where $\tilde{p}_{I\rightarrow i}(w_{i})\propto{m}_{I\rightarrow i}(w_{i})$.

Our proposed approach is to regard Eq. 5 as expectation value of
$f_{I}(\underline{w}_{I})$ with respect to the joint probability distribution
$p({w}_{\partial{I}\setminus{i}})$. The value of the message can be estimated
using a Monte Carlo (MC) approach, _i.e._ , drawing samples from the joint
distribution:

$\tilde{p}_{I\rightarrow i}(w_{i};p_{\partial I\setminus i\rightarrow
I})=\mathbb{E}_{p({w}_{\partial{I}\setminus{i}})}\left[f_{I}(\underline{w}_{I})\right]\approx\frac{1}{L}\sum_{m=1}^{L}f_{I}({\underline{w}_{I,m}}),$
(6)

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}\in\\{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 $i$th 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.

Input: Normalized messages $\\{{p}_{i\rightarrow I}(w_{i})\\}_{\\{i\in\partial
I\\}}$ incident on factor $I$

Output: $\\{{p}_{I\rightarrow i}(w_{i}=1)\\}_{\\{i\in\partial I\\}}$,
$\tilde{p}_{I\rightarrow i}(w_{i}=1)$ and $\tilde{p}_{I\rightarrow
i}(w_{i}=0)$ from factor $I$

1 for  _$i\in\partial I$_ do

2 for  _$m\in\\{1,\dots,L\\}$_ do

3 $\underline{w}^{1}_{I,m}\sim\\{\prod_{j\in\partial{I}}{p}_{j\rightarrow
I}(w_{j})\\}\quad$

4 ${w}^{1}_{i,m}\leftarrow 1\quad$

5 $\underline{w}^{0}_{I,m}\sim\\{\prod_{j\in\partial{I}}{p}_{j\rightarrow
I}(w_{j})\\}\quad$

6 ${w}^{0}_{i,m}\leftarrow 0\quad$

7 end for

8

9 $\tilde{p}_{I\rightarrow
i}(w_{i}=1)\leftarrow\frac{1}{L}\sum_{m=1}^{l}f_{I}({\underline{w}^{1}_{I,m}})$

10 $\tilde{p}_{I\rightarrow
i}(w_{i}=0)\leftarrow\frac{1}{L}\sum_{m=1}^{l}f_{I}({\underline{w}^{0}_{I,m}})$

11 ${p}_{I\rightarrow i}(w_{i}=1)\leftarrow\frac{\tilde{p}_{I\rightarrow
i}(w_{i}=1)}{\tilde{p}_{I\rightarrow i}(w_{i}=0)+\tilde{p}_{I\rightarrow
i}(w_{i}=1)}$

12 end for

Algorithm 1 Stochastic BP factor to variable update subroutine for binary
variables: StochasticBPF2Var$(\\{{p}_{i\rightarrow
I}(w_{i})\\}_{\\{i\in\partial I\\}})$.

###  3.2 Stochastic SP (S3P and S4P)

The integrals in the SP message updates are converted into sums by keeping
histograms of BP messages. With $K$ bins, the sums involve $K^{|\partial
I|-1}$ and $K^{|\partial 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 $\mathcal{O}(K^{2}|\partial 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 $\approx|\partial i|K$ in order to avoid aliasing. Since
$|\partial 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 defining
$\tilde{S}_{j\rightarrow I}({p}_{j\rightarrow I})\propto{{M}_{j\rightarrow
I}({p}_{j\rightarrow I})}$ with ${S}_{j\rightarrow I}({p}_{j\rightarrow
I})\triangleq\frac{\tilde{S}_{j\rightarrow I}({p}_{j\rightarrow
I})}{\sum_{{p}_{j\rightarrow I}}\tilde{S}_{j\rightarrow I}({p}_{j\rightarrow
I})}$ , we can regard the factor to variable update as an expectation under
the joint distribution $S({p}_{\partial I\setminus
i})\triangleq\prod_{j\in\partial I\setminus i}{S}_{j\rightarrow
I}({p}_{j\rightarrow I})$. Note that each ${S}_{j\rightarrow
I}({p}_{j\rightarrow I})$ is a categorical distribution (with cardinality
$K$). The resulting factor to variable message is:

$\tilde{S}_{I\rightarrow i}({p}_{I\rightarrow
i}=p_{k})=\mathbb{E}_{S({p}_{\partial I\setminus
i})}\Bigg{[}\mathbb{I}\left[p_{k}={p}_{I\rightarrow i}(.;p_{\partial
I\setminus i\rightarrow I})\right]{p}_{I\rightarrow i}(.;p_{\partial
I\setminus i\rightarrow I})(\emptyset)\Bigg{]},$ (7)

where ${p}_{k}\in[0,1]$ is one of the possible $K$ values that
${p}_{I\rightarrow i}$ can take. ${p}_{I\rightarrow i}(.;p_{\partial
I\setminus i\rightarrow I})(\emptyset)$ 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\rightarrow
I}({p}_{{i\rightarrow I}})\\}_{\\{i\in\partial 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
$\texttt{bin}(p_{k})\triangleq 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 S4P (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 S3P
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 S4P the Stochastic Nonlocal Message Passing (SNMP) model. SNMP starts by
running the faster BP first, and if all the constraints are not satisfied, S4P
is called.

Input: Normalized messages $\\{{S}_{i\rightarrow I}({p}_{{i\rightarrow
I}})\\}_{\\{i\in\partial I\\}}$ incident on factor $I$

Output: Normalized outgoing messages $\\{{S}_{I\rightarrow
i}({p}_{{i\rightarrow I}})\\}_{\\{i\in\partial I\\}}$ from factor $I$

1 for  _$m\in\\{1,\dots L_{\text{sp}}\\}$_ do

2 ${p}_{\partial I\rightarrow I}\sim\prod_{j\in\partial I}{S}_{j\rightarrow
I}({p}_{j\rightarrow I})$

3 ${p}_{I\rightarrow i}(w_{\partial I}=1),{\tilde{p}}_{I\rightarrow
i}(w_{\partial I}=1),{\tilde{p}}_{I\rightarrow i}(w_{\partial I}=0)$
$\leftarrow$ StochasticBPF2Var$({p}_{\partial I\rightarrow I})$

4 for  _$i\in\partial I$_ do

5

6 ${p}_{I\rightarrow i}(\emptyset)\leftarrow{\tilde{p}}_{I\rightarrow
i}(w_{i}=0)+{\tilde{p}}_{I\rightarrow i}(w_{i}=1)$

7 ${k}\leftarrow\texttt{bin}({p}_{I\rightarrow i}(w_{i}=1))$

8 $\tilde{S}_{I\rightarrow i}(p_{k})\leftarrow\tilde{S}_{I\rightarrow
i}(p_{k})+{p}_{I\rightarrow i}(\emptyset)$

9 end for

10

11 end for

12for  _$i\in\partial I$_ do

13 for  _$k\in\\{1,...,K\\}$_ do

14 ${S}_{I\rightarrow i}({p}_{k})\leftarrow\frac{\tilde{S}_{I\rightarrow
i}({p}_{k})}{\sum_{k=1}^{K}\tilde{S}_{I\rightarrow i}({p}_{k})}$

15 end for

16

17 end for

Algorithm 2 SP Stochastic factor to variable subroutine.

###  3.3 Problem set-up

Consider a training set of binary vectors $\underline{x}$ and labels,
$\underline{y}$, comprising of $M$ examples. We assume binary labels,
$y_{I}\in\\{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 $\underline{w}$ that predicts label $g_{\underline{w}}(x_{I})$ given
input $x_{I}$. $g_{\underline{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\triangleq|\underline{w}|$
variables: $F(\underline{w})=\prod_{I=1}^{M}f_{I}(\underline{w})$ and
$f_{I}(\underline{w})\triangleq\mathbb{I}\left[g_{\underline{w}}(x_{I})=y_{I}\right]$.
We have considered the general case where each factor $f_{I}$ depends on all
the variables, _i.e._ , $\underline{w}_{I}=\underline{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(\underline{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
$F(\underline{w})\propto\exp(-E(\underline{w}))$, where the energy is defined
as $E\triangleq-\sum_{I}E_{I}(\underline{w}))$ with
$E_{I}(\underline{w})\triangleq-\beta(1-f_{I}(\underline{w}))$. Here, $\beta$
is the inverse temperature for a Boltzmann distributions in the language of
statistical mechanics (assuming Boltzmann constant is unity). A high value of
$\beta$ penalizes the unsatisfied factors more severely.

##  4 Experiments

We first demonstrate in Fig. 1 that SBP and S4P yield comparable results to
the BP and S3P, 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.

Figure 1: (left) Comparing the performance of stochastic and deterministic BP.
The red dashed line denotes a line with unit slope, and the size of each
marker denotes the number of samples ($L_{BP}$) used in the MC estimate for
BP. We find the algorithms are comparable even with very small $L_{BP}$.
(right) Performance of S4P against S3P: for larger values of samples S4P is
able to reach higher accuracies, due to accumulative binning errors that could
arise in S3P (see Appendix. A for details).

###  4.1 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 $\underline{x}$ and $\underline{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
$\alpha\triangleq\frac{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 $\alpha$), many
solutions exist. For an under-parameterized system at large values of
$\alpha$, 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 $\alpha$ 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 $\alpha$ 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, $\underline{y}$, have a
functional relationship with $\underline{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.

Figure 2: Comparing the performance of the proposed models and GD on several
datasets under different settings (see text for detail).

####  4.1.1 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.

Figure 3: Visualizing the evolution of the stochastic SP marginals, for
different sample sizes ($L_{sp}$). 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
S4P can find better configurations which are inaccessible to BP. Best viewed
in colour.

###  4.2 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). Fig. 5 compares
the results of running SNMP with different values of $\beta$ with the SGD
algorithm. The best performance is achieved at higher values of $\beta$, but
not at $\beta=\infty$. The results are comparable with SGD.

Figure 4: Comparing performance of the different solvers on the MNIST dataset.
The two leftmost plots correspond respectively to the training and test
performance of a linear model, whereas the latter two are obtained from a two-
layer neural network with sign-function activations. 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 $\beta$ leads to a
higher accuracy score, though noticeably $\beta_{\infty}$ 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.

##  5 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.

## References

  * Achterhold et al. [2018] Jan Achterhold, Jan M. Köhler, Anke Schmeink, and Tim Genewein.  Variational network quantization.  In _ICLR_ , 2018.
  * Alpaydin and Alimoglu [1998] E. Alpaydin and Fevzi. Alimoglu.  Pen-based recognition of handwritten digits data set.  https://archive.ics.uci.edu/ml/datasets/Pen-Based+Recognition+of+Handwritten+Digits, 1998\.  Accessed: 2019-09-22.
  * Bishop [2006] Christopher M Bishop.  _Pattern recognition and machine learning_.  springer, 2006.
  * Braunstein and Zecchina [2006] Alfredo Braunstein and Riccardo Zecchina.  Learning by message passing in networks of discrete synapses.  _Physical review letters_ , 96(3):030201, 2006\.
  * Braunstein et al. [2005] Alfredo Braunstein, Marc Mézard, and Riccardo Zecchina.  Survey propagation: An algorithm for satisfiability.  _Random Structures & Algorithms_, 27(2):201–226, 2005.
  * Chieu and Lee [2009] Hai Leong Chieu and Wee Sun Lee.  Relaxed survey propagation for the weighted maximum satisfiability problem.  _J. Artif. Intell. Res._ , 36:229–266, 2009.
  * Courbariaux and Bengio [2016] Matthieu Courbariaux and Yoshua Bengio.  Binarynet: Training deep neural networks with weights and activations constrained to +1 or -1.  _CoRR_ , abs/1602.02830, 2016.  URL http://arxiv.org/abs/1602.02830.
  * Heskes [2002] Tom Heskes.  Stable fixed points of loopy belief propagation are local minima of the bethe free energy.  In _NIPS_ , 2002.
  * Koller et al. [2009] Daphne Koller, Nir Friedman, and Francis Bach.  _Probabilistic graphical models: principles and techniques_.  MIT press, 2009.
  * Kschischang et al. [2001] Frank R Kschischang, Brendan J Frey, Hans-Andrea Loeliger, et al.  Factor graphs and the sum-product algorithm.  _IEEE Transactions on information theory_ , 47(2):498–519, 2001.
  * Louizos et al. [2017] Christos Louizos, Karen Ullrich, and Max Welling.  Bayesian compression for deep learning.  In _NIPS_ , 2017.
  * Mezard and Montanari [2009] Marc Mezard and Andrea Montanari.  _Information, physics, and computation_.  Oxford University Press, 2009.
  * Murphy et al. [2013] Kevin P. Murphy, Yair Weiss, and Michael I. Jordan.  Loopy belief propagation for approximate inference: An empirical study.  _CoRR_ , abs/1301.6725, 2013.  URL http://arxiv.org/abs/1301.6725.
  * Noorshams and Wainwright [2013] Nima Noorshams and Martin J. Wainwright.  Stochastic belief propagation: A low-complexity alternative to the sum-product algorithm.  _IEEE Transactions on Information Theory_ , 59(4):1981–2000, Apr 2013.  ISSN 1557-9654.  doi: 10.1109/tit.2012.2231464.  URL http://dx.doi.org/10.1109/TIT.2012.2231464.
  * Paszke et al. [2017] Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang, Zachary DeVito, Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer.  Automatic differentiation in pytorch.  2017\.
  * Peters and Welling [2018] Jorn WT Peters and Max Welling.  Probabilistic binary neural networks.  _arXiv preprint arXiv:1809.03368_ , 2018.
  * Ravanbakhsh and Greiner [2014] Siamak Ravanbakhsh and Russell Greiner.  Revisiting algebra and complexity of inference in graphical models.  _arXiv preprint arXiv:1409.7410_ , 2014.
  * Shayer et al. [2018] Oran Shayer, Dan Levi, and Ethan Fetaya.  Learning discrete weights using the local reparameterization trick.  _CoRR_ , abs/1710.07739, 2018.
  * Srinivasa et al. [2016] Christopher Srinivasa, Siamak Ravanbakhsh, and Brendan Frey.  Survey propagation beyond constraint satisfaction problems.  In Arthur Gretton and Christian C. Robert, editors, _Proceedings of the 19th International Conference on Artificial Intelligence and Statistics_ , volume 51 of _Proceedings of Machine Learning Research_ , pages 286–295, Cadiz, Spain, 09–11 May 2016. PMLR.  URL http://proceedings.mlr.press/v51/srinivasa16.html.
  * Srinivasa et al. [2017] Christopher Srinivasa, Inmar Givoni, Siamak Ravanbahksh, and Brendan J. Frey.  Min-max propagation.  2017\.  URL https://papers.nips.cc/paper/7140-min-max-propagation.
  * Tarlow et al. [2010] Daniel Tarlow, Inmar Givoni, and Richard Zemel.  Hop-map: Efficient message passing with high order potentials.  In Yee Whye Teh and Mike Titterington, editors, _Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics_ , volume 9 of _Proceedings of Machine Learning Research_ , pages 812–819, Chia Laguna Resort, Sardinia, Italy, 13–15 May 2010. PMLR.  URL http://proceedings.mlr.press/v9/tarlow10a.html.
  * Ullrich et al. [2017] Karen Ullrich, Edward Meeds, and Max Welling.  Soft weight-sharing for neural network compression.  _CoRR_ , abs/1702.04008, 2017.
  * Yedidia et al. [2005] Jonathan S Yedidia, William T Freeman, and Yair Weiss.  Constructing free-energy approximations and generalized belief propagation algorithms.  _IEEE Transactions on information theory_ , 51(7):2282–2312, 2005.
  * Zhang et al. [2017] Chiyuan Zhang, Samy Bengio, Moritz Hardt, Benjamin Recht, and Oriol Vinyals.  Understanding deep learning requires rethinking generalization.  _CoRR_ , abs/1611.03530, 2017.

Supplementary material

## Supplementary Material for Nonlocal optimization of binary neural networks

##  Appendix A Consistency between stochastic and exact models

We gauge the accuracy of SBP and S4P against BP and S3P on a large number of
problems and range of parameters for which BP and S3P are tractable. More
specifically, we fix $M=2^{N}$, and sample all combinations from a grid formed
by:

  * •

$N\in\\{4,6,8,10\\}$

  * •

$\gamma\in\\{0.2,0.5,0.8\\}$

  * •

$L_{BP}\in\\{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.

In a similar way, to compare S3P and S4P, we fix $N=20$ and sample all
combinations from a grid of:

  * •

$K\in\\{31,51,101,201\\}$

  * •

$\gamma\in\\{0.2,0.5,0.8\\}$

  * •

$L_{BP}\in\\{5,10,20,50,100\\}$

  * •

$L_{SP}\in\\{5,10,20,50,100\\}$

  * •

$M\in\\{8,12,16,32\\}$

A comparison between these accuracies is also present on the right side of
Fig. 1. Two important points are noted: first, higher numbers of MC samples
lead to better S4P performance. Second, quite interestingly, S4P usually
outperforms S3P, due to the binning errors that accumulate if all the the
exponential combinations of bins are considered.

##  Appendix 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:

$\tilde{S}_{i\rightarrow I}({p}_{i\rightarrow
I}=p_{k})=\mathbb{E}_{S({p}_{\partial i\setminus
I})}\Bigg{[}\mathbb{I}\left[p_{k}={p}_{i\rightarrow I}(.;p_{\partial
i\setminus I\rightarrow i})\right]{p}_{i\rightarrow I}(.;p_{\partial
i\setminus I\rightarrow i})(\emptyset)\Bigg{]},$ (8)

where ${p}_{k}\in[0,1]$ is one of the possible $K$ values that
${p}_{i\rightarrow I}$ can take. ${p}_{I\rightarrow i}(.;p_{\partial
I\setminus i\rightarrow I})(\emptyset)$ 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\rightarrow
I}({p}_{{i\rightarrow I}})\\}_{\\{i\in\partial 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
$\texttt{bin}(p_{k})\triangleq k$. The message corresponding to this bin is
subsequently updated (line 7). Finally, the normalized messages are calculated
(lines 10-14).

Input: Normalized messages $\\{{S}_{I\rightarrow i}({p}_{{I\rightarrow
i}})\\}_{\\{I\in\partial i\\}}$ incident on variable $i$

Output: Normalized outgoing messages $\\{{S}_{i\rightarrow
I}({p}_{{i\rightarrow I}})\\}_{\\{i\in\partial I\\}}$ from variable $i$

1 for  _$m\in\\{1,\dots L_{\text{sp}}\\}$_ do

2 ${p}_{\partial i\rightarrow i}\sim\prod_{J\in\partial i}{S}_{J\rightarrow
i}({p}_{J\rightarrow i})$

3 ${p}_{i\rightarrow I}(w_{i}=1),{\tilde{p}}_{i\rightarrow
I}(w_{i}=1),{\tilde{p}}_{i\rightarrow I}(w_{i}=0)$ $\leftarrow$
BPVar2F$({p}_{\partial i\rightarrow i})$

4 for  _$I\in\partial i$_ do

5

6 ${p}_{i\rightarrow I}(\emptyset)\leftarrow{\tilde{p}}_{i\rightarrow
I}(w_{i}=0)+{\tilde{p}}_{i\rightarrow I}(w_{i}=1)$

7 ${k}\leftarrow\texttt{bin}({p}_{i\rightarrow I}(w_{i}=1))$

8 $\tilde{S}_{i\rightarrow I}(p_{k})\leftarrow\tilde{S}_{i\rightarrow
I}(p_{k})+{p}_{i\rightarrow I}(\emptyset)$

9 end for

10

11 end for

12for  _$I\in\partial i$_ do

13 for  _$k\in\\{1,...,K\\}$_ do

14 ${S}_{i\rightarrow I}({p}_{k})\leftarrow\frac{\tilde{S}_{i\rightarrow
I}({p}_{k})}{\sum_{k=1}^{K}\tilde{S}_{i\rightarrow I}({p}_{k})}$

15 end for

16

17 end for

Algorithm 3 SP Stochastic variable to factor subroutine.

##  Appendix 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_{BP}$: the number of samples used for SBP. As shown on the left of Fig. 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 S3P and S4P. 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.

  * •

$\gamma$: 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 $\mu_{t}$, then with damping
factor $\gamma$ after computing updates $\mu$, the new message becomes:
$\mu_{t+1}=(1-\gamma)\mu_{t}+\mu$. When $L_{SP}$ is high, $\gamma$ can afford
to be higher, and converges more quickly. For SBP we typically set
$\gamma=0.2$, while for S4P we use $\gamma=0.8$.

  * •

$\beta$: the inverse temperature for Max-Sat. Higher values of $\beta$ have a
positive impact on the variance encountered during training, especially when
combined with tuning $\gamma$. Fig. 5 provides a good illustration of the
impact of $\beta$ 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$.

For both SBP and S4P, we set $L_{BP}=5$. For S4P, we use $L_{SP}=100$ for all
experiments, with $K=201$. In all our experiments, we keep the number of
variables fixed and adjust the number of training examples.

For the glass dataset, we set $N=10$ and let $M$ span $[5...50]$ jumping in
units of 5. We use the same span of $M$ for stained glass, but set $N=31$. For
the digits dataset, the linear model requires $N=64$ to match the input
dimensions. For $M$, we span the range $[5...245]$ jumping in units of 10.

For Fig. 3, we set $N=10$, $M=20$, and ran for 20 iterations with
$\gamma=0.5$. $K$ was set to $31$ for ease of visualization.

###  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\times 14=196$. The
linear model thus consists of 196 parameters, while the MLP, with a single
hidden layer of size 2, has $196\times 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 $\gamma=0.5$ for the linear model and the MLP. For the convolutional
network, we adapt $\gamma$ according to the present accuracy of the model. At
the beginning of training, when accuracy is low, $\gamma$ should be high to
explore the space quickly. As the model becomes more accurate, however,
$\gamma$ can be reduced to limit the variance. As such, every 20 iterations we
compute the accuracy, and set $\gamma=(1-\textrm{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 $\beta$
in the range $\\{0.1,1,10,100,\infty\\}$, where $\infty$ 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 S3P and S4P 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 $[B\times C_{out},C_{in},K_{W},K_{H}]$, where $B=M\times
L_{SP}\times L_{BP}\times 2N$. The inputs to the layer are reshaped to
$[1,B\times 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.