Year: 2021

Abstract:
Hopfield networks (HNs) and Restricted Boltzmann Machines (RBMs) are two important models at the interface of statistical physics, machine learning, and neuroscience. Recently, there has been interest in the relationship between HNs and RBMs, due to their similarity under the statistical mechanics formalism. An exact mapping between HNs and RBMs has been previously noted for the special case of orthogonal (“uncorrelated”) encoded patterns. We present here an exact mapping in the case of correlated pattern HNs, which are more broadly applicable to existing datasets. Specifically, we show that any HN with N binary variables and p < N arbitrary binary patterns can be transformed into an RBM with N binary visible variables and p gaussian hidden variables. We outline the conditions under which the reverse mapping exists, and conduct experiments on the MNIST dataset which suggest the mapping provides a useful initialization to the RBM weights. We discuss extensions, the potential importance of this correspondence for the training of RBMs, and for understanding the performance of deep architectures which utilize RBMs.

heading: 1 INTRODUCTION
text: Hopfield networks (HNs) (Hopfield, 1982; Amit, 1989) are a classical neural network architecture that can store prescribed patterns as fixed-point attractors of a dynamical system. In their standard formulation with binary valued units, HNs can be regarded as spin glasses with pairwise interactions Jij that are fully determined by the patterns to be encoded. HNs have been extensively studied in the statistical mechanics literature (e.g. (Kanter & Sompolinsky, 1987; Amit et al., 1985)), where they can be seen as an interpolation between the ferromagnetic Ising model (p = 1 pattern) and the Sherrington-Kirkpatrick spin glass model (many random patterns) (Kirkpatrick & Sherrington, 1978; Barra & Guerra, 2008). By encoding patterns as dynamical attractors which are robust to perturbations, HNs provide an elegant solution to pattern recognition and classification tasks. They are considered the prototypical attractor neural network, and are the historical precursor to modern recurrent neural networks.
Concurrently, spin glasses have been used extensively in the historical machine learning literature where they comprise a sub-class of “Boltzmann machines” (BMs) (Ackley et al., 1985). Given a collection of data samples drawn from a data distribution, one is generally interested in “training” a BM by tuning its weights Jij such that its equilibrium distribution can reproduce the data distribution as closely as possible (Hinton, 2012). The resulting optimization problem is dramatically simplified when the network has a two-layer structure where each layer has no self-interactions, so that there are only inter-layer connections (Hinton, 2012) (see Fig. 1). This architecture is known as a Restricted Boltzmann Machine (RBM), and the two layers are sometimes called the visible layer and the hidden layer. The visible layer characteristics (dimension, type of units) are determined by the training data, whereas the hidden layer can have binary or continuous units and the dimension is chosen somewhat arbitrarily. In addition to generative modelling, RBMs and their multi-layer extensions have been used for a variety of learning tasks, such as classification, feature extraction, and dimension reduction (e.g. Salakhutdinov et al. (2007); Hinton & Salakhutdinov (2006)).
There has been extensive interest in the relationship between HNs and RBMs, as both are built on the Ising model formalism and fulfill similar roles, with the aim of better understanding RBM behaviour and potentially improving performance. Various results in this area have been recently reviewed (Marullo & Agliari, 2021). In particular, an exact mapping between HNs and RBMs has been previously noted for the special case of uncorrelated (orthogonal) patterns (Barra et al., 2012). Several related models have since been studied (Agliari et al., 2013; Mézard, 2017), which partially relax the uncorrelated pattern constraint. However, the patterns observed in most real datasets exhibit significant correlations, precluding the use of these approaches.
In this paper, we demonstrate exact correspondence between HNs and RBMs in the case of correlated pattern HNs. Specifically, we show that any HN with N binary units and p < N arbitrary (i.e. non-orthogonal) binary patterns encoded via the projection rule (Kanter & Sompolinsky, 1987; Personnaz et al., 1986), can be transformed into an RBM with N binary and p gaussian variables. We then characterize when the reverse map from RBMs to HNs can be made. We consider a practical example using the mapping, and discuss the potential importance of this correspondence for the training and interpretability of RBMs.

heading: 2 RESULTS
text: We first introduce the classical solution to the problem of encodingN -dimensional binary {−1,+1} vectors {ξµ}pµ=1, termed “patterns”, as global minima of a pairwise spin glass H(s) = − 12s
TJs. This is often framed as a pattern retrieval problem, where the goal is to specify or learn Jij such that an energy-decreasing update rule for H(s) converges to the patterns (i.e. they are stable fixed points). Consider the N × p matrix ξ with the p patterns as its columns. Then the classical prescription known as the projection rule (or pseudo-inverse rule) (Kanter & Sompolinsky, 1987; Personnaz et al., 1986), J = ξ(ξT ξ)−1ξT , guarantees that the p patterns will be global minima of H(s). This resulting spin model is commonly called a (projection) Hopfield network, and has the Hamiltonian
H(s) = −1 2 sT ξ(ξT ξ)−1ξTs. (1)
Note that ξT ξ invertibility is guaranteed as long as the patterns are linearly independent (we therefore require p ≤ N ). Also note that in the special (rare) case of orthogonal patterns ξµ · ξν = Nδµν (also called “uncorrelated”), studied in the previous work (Barra et al., 2012), one has ξT ξ = NI and so the pseudo-inverse interactions reduce to the well-known Hebbian form J = 1N ξξ
T (the properties of which are studied extensively in Amit et al. (1985)). Additional details on the projection HN Eq. (1) are provided in Appendix A. To make progress in analyzing Eq. (1), we first consider a transformation of ξ which eliminates the inverse factor.

heading: 2.1 MAPPING A HOPFIELD NETWORK TO A RESTRICTED BOLTZMANN MACHINE
text: In order to obtain a more useful representation of the quadratic form Eq. (1) (for our purposes), we utilize the QR-decomposition (Schott & Stewart, 1999) of ξ to “orthogonalize” the patterns,
ξ = QR, (2)
with Q ∈ RN×p,R ∈ Rp×p. The columns of Q are the orthogonalized patterns, and form an orthonormal basis (of non-binary vectors) for the p-dimensional subspace spanned by the binary patterns. R is upper triangular, and if its diagonals are held positive then Q and R are both unique (Schott & Stewart, 1999). Note both the order and sign of the columns of ξ are irrelevant for HN pattern recall, so there are n = 2p · p! possibleQ,R pairs. Fixing a pattern ordering, we can use the orthogonality ofQ to re-write the interaction matrix as
J = ξ(ξT ξ)−1ξT = QR(RTR)−1RTQT = QQT (3)
(the last equality follows from (RTR)−1 = R−1(RT )−1). Eq. (3) resembles the simple Hebbian rule but with non-binary orthogonal patterns. Defining q ≡ QTs in analogy to the classical pattern overlap parameterm ≡ 1N ξ Ts (Amit et al., 1985), we have
H(s) = −1 2 sTQQTs = −1 2 q(s) · q(s). (4)
Using a Gaussian integral as in Amit et al. (1985); Barra et al. (2012); Mézard (2017) to transform (exactly) the partition function Z ≡ ∑ {s} e −βH(s) of Eq. (1), we get
Z = ∑ {s} e 1 2 (βq) T (β−1I)(βq)
= ∑ {s} ∫ e− β 2 ∑ µ λ 2 µ+β ∑ µ λµ ∑ iQiµsi ∏ µ dλµ√ 2π/β . (5)
The second line can be seen as the partition function of an expanded Hamiltonian for the N (binary) original variables {si} and the p (continuous) auxiliary variables {λµ}, i.e.
HRBM({si}, {λµ}) = 1
2 ∑ µ λ2µ − ∑ µ ∑ i Qiµsiλµ. (6)
Note that this is the Hamiltonian of a binary-continuous RBM with inter-layer weights Qiµ. The original HN is therefore equivalent to an RBM described by Eq. (6) (depicted in Fig. 1). As mentioned above, there are many RBMs which correspond to the same HN due to the combinatorics of choosing Q. In fact, instead of QR factorization one can use any decomposition which satisfies J = UUT , with orthogonal U ∈ RN×p (see Appendix B), in which case U acts as the RBM weights. Also note the inclusion of an applied field term − ∑ i bisi in Eq. (1) trivially carries
through the procedure, i.e. H̃RBM ({si}, {λµ}) = 12 ∑ µ λ 2 µ − ∑ i bisi − ∑ µ ∑ iQiµsiλµ.
Instead of working with the joint form Eq. (6), one could take a different direction from Eq. (5) and sum out the original variables {si}, i.e.
Z = ∫ e− β 2 ∑ µ λ 2 µ2N ∏ i cosh ( β ∑ µ Qiµλµ )∏ µ dλµ√ 2π/β . (7)
This continuous, p-dimensional representation is useful for numerical estimation of Z (Section 3.1). We may write Eq. (7) as Z = ∫ e−F0(λ)dλµ, where
F0({λµ}) = 1
2 ∑ µ λ2µ − 1 β ∑ i ln cosh ( β ∑ µ Qiµλµ ) . (8)
Eq. (8) is an approximate Lyapunov function for the mean dynamics of {λµ}; ∇λF0 describes the effective behaviour of the stochastic dynamics of the N binary variables {si} at temperature β−1.

heading: 2.2 COMMENTS ON THE REVERSE MAPPING
text: With the mapping from HNs (with correlated patterns) to RBMs established, we now consider the reverse direction. Consider a binary-continuous RBM with inter-layer weights Wiµ which couple a
visible layer of N binary variables {si} to a hidden layer of p continuous variables {λµ},
H(s,λ) = 1
2 ∑ µ λ2µ − ∑ i bisi − ∑ µ ∑ i Wiµsiλµ. (9)
Here we use W instead of Q for the RBM weights to emphasize that the RBM is not necessarily an HN. First, following Mehta et al. (2019), we transform the RBM to a BM with binary states by integrating out the hidden variables. The corresponding Hamiltonian for the visible units alone is (see Appendix D.1 for details),
H̃(s) = − ∑ i bisi − 1 2 ∑ i ∑ j ∑ µ WiµWjµsisj , (10)
a pairwise Ising model with a particular coupling structure Jij = ∑ µWiµWjµ, which in vector form is J =
∑ µ wµw T µ =WW T , (11)
where {wµ} are the p columns ofW . In general, this Ising model Eq. (10) produced by integrating out the hidden variables need not have Hopfield structure (discussed below). However, it automatically does (as noted in Barra et al. (2012)), in the very special case whereWiµ ∈ {−1,+1}. In that case, the binary patterns are simply {wµ}, so that Eq. (11) represents a Hopfield network with the Hebbian prescription. This situation is likely rare and may only arise as a by-product of constrained training; for a generically trained RBM the weights will not be binary. It is therefore interesting to clarify when and how real-valued RBM interactionsW can be associated with HNs.
Approximate binary representation of W : In Section 2.1, we orthogonalized the binary matrix ξ via the QR decomposition ξ = QR, where Q is an orthogonal (but non-binary) matrix, which allowed us to map a projection HN (defined by its patterns ξ, Eq. (1)) to an RBM (defined by its inter-layer weightsQ, Eq. (6)).
Here we consider the reverse map. Given a trained RBM with weights W ∈ RN×p, we look for an invertible transformation X ∈ Rp×p which binarizes W . We make the mild assumption that W is rank p. If we find such an X , then B =WX will be the Hopfield pattern matrix (analogous to ξ), with Biµ ∈ {−1,+1}. This is a non-trivial problem, and an exact solution is not guaranteed. As a first step to study the problem, we relax it to that of finding a matrix X ∈ GLp(R) (i.e. invertible, p × p, real) which minimizes the binarization error
argmin X∈GLp(R)
||WX − sgn(WX)||F . (12)
We denote the approximately binary transformation ofW via a particular solutionX by
Bp =WX. (13)
We also define the associated error matrixE ≡ Bp− sgn(Bp). We stress thatBp is non-binary and approximatesB ≡ sgn(Bp), the columns of which will be HN patterns under certain conditions on E. We provide an initial characterization and example in Appendix D.

heading: 3 EXPERIMENTS ON MNIST DATASET
text: Next we investigate whether the Hopfield-RBM correspondence can provide an advantage for training binary-gaussian RBMs. We consider the popular MNIST dataset of handwritten digits (LeCun et al., 1998) which consists of 28× 28 images of handwritten images, with greyscale pixel values 0 to 255. We treat the sample images as N ≡ 784 dimensional binary vectors of {−1,+1} by setting all non-zero values to +1. The dataset includes M ≡ 60, 000 training images and 10, 000 testing images, as well as their class labels µ ∈ {0, . . . , 9}.

heading: 3.1 GENERATIVE OBJECTIVE
text: The primary task for generative models such as RBMs is to reproduce a data distribution. Given a data distribution pdata, the generative objective is to train a model (here an RBM defined by its parameters θ), such that the model distribution pθ is as close to pdata as possible. This is often quantified by the Kullback-Leibler (KL) divergence DKL(pdata‖pθ) = ∑ s pdata(s) ln ( pdata(s) pθ(s) ) . One generally does not have access to the actual data distribution, instead there is usually a representative training set S = {sa}Ma=1 sampled from it. As the data distribution is constant with respect to θ, the generative objective is equivalent to maximizing L(θ) = 1M ∑ a ln pθ(sa).

heading: 3.1.1 HOPFIELD RBM SPECIFICATION
text: With labelled classes of training data, specification of an RBM via a one-shot Hopfield rule (“Hopfield RBM”) is straightforward. In the simplest approach, we define p = 10 representative patterns via the (binarized) class means
ξµ ≡ sgn  1 |Sµ| ∑ s∈Sµ s  . (14) where µ ∈ {0, . . . , 9} and Sµ is the set of sample images for class µ. These patterns comprise the columns of theN×p pattern matrix ξ, which is then orthogonalized as in Eq. (2) to obtain the RBM weightsW which couple N binary visible units to p gaussian hidden units.
We also consider refining this approach by considering sub-classes within each class, representing, for example, the different ways one might draw a “7”. As a proof of principle, we split each digit class into k sub-patterns using hierarchical clustering. We found good results with Agglomerative clustering using Ward linkage and Euclidean distance (see Murtagh & Contreras (2012) for an overview of this and related methods). In this way, we can define a hierarchy of Hopfield RBMs. At one end, k = 1, we have our simplest RBM which has p = 10 hidden units and encodes 10 patterns (using Eq. (14)), one for each digit class. At the other end, 10k/N → 1, we can specify increasingly refined RBMs that encode k sub-patterns for each of the 10 digit classes, for a total of p = 10k patterns and hidden units. This approach has an additional cost of identifying the sub-classes, but is still typically faster than training the RBM weights directly (discussed below).
The generative performance as a function of k and β is shown in Fig. 2, and increases monotonically with k in the range plotted. If β is too high (very low temperature) the free energy basins will be very deep directly at the patterns, and so the model distribution will not capture the diversity of images from the data. If β is too low (high temperature), there is a “melting transition” where the original pattern basins disappear entirely, and the data will therefore be poorly modelled. Taking α = p/N ∼ 0.1 (roughly k = 8), Fig. 1 of Kanter & Sompolinsky (1987) predicts βm ≈ 1.5 for the theoretical melting transition for the pattern basins. Interestingly, this is quite close to our observed peak near β = 2. Note also as k is increased, the generative performance is sustained at lower temperatures.
In situations where one already has access to the class labels, this approach to obtain RBM weights is very fast. The class averaging has negligible computational cost O(MN) for the whole training
set (M samples), and the QR decomposition has a modest complexity ofO(Np2) (Schott & Stewart, 1999). Conventional RBM training, discussed below, requires significantly more computation.

heading: 3.1.2 CONVENTIONAL RBM TRAINING
text: RBM training is performed through gradient ascent on the log-likelihood of the data, L(θ) = 1 M ∑ a ln pθ(sa) (equivalent here to minimizing KL divergence, as mentioned above). We are focused here on the weights W in order to compare to the Hopfield RBM weights, and so we neglect the biases on both layers. As is common (Hinton, 2012), we approximate the total gradient by splitting the training dataset into “mini-batches”, denoted B. The resulting gradient ascent rule for the weights is (see Appendix E)
W t+1iµ =W t iµ + η 〈∑ j Wjµs a i s a j 〉 a∈B − 〈siλµ〉model  , (15) where 〈siλµ〉model ≡ Z−1 ∑ {s} ∫ siλµe −βH(s,λ)∏ µ dλµ is an average over the model distribution.
The first bracketed term of Eq. (15) is simple to calculate at each iteration of the weights. The second term, however, is intractable as it requires one to calculate the partition function Z. We instead approximate it using contrastive divergence (CD-K) (Carreira-Perpinan & Hinton, 2005; Hinton, 2012). See Appendix E for details. Each full step of RBM weight updates involves O(KBNp) operations (Melchior et al., 2017). Training generally involves many mini-batch iterations, such that the entire dataset is iterated over (one epoch) many times. In our experiments we train for 50 epochs with mini-batches of size 100 (3 · 105 weight updates), so the overall training time can be extensive compared to the one-shot Hopfield approach presented above. For further details on RBM training see e.g. Hinton (2012); Melchior et al. (2017).
In Fig. 3, we give an example of the Hopfield RBM weights (for k = 1), as well as how they evolve during conventional RBM training. Note Fig. 3(a), (b) appear qualitatively similar, suggesting that the proposed initializationQ from Eqs. (2), (14) may be near a local optimum of the objective.
In Fig. 4(a), (b), we compare conventional RBM training on four different weight initializations: (i) randomWiµ ∼ N (0, 0.01) (purple), commonly used in the literature; (ii) our proposed weights from the projection rule Hopfield mapping for correlated patterns (blue); (iii) the “Hebbian” Hopfield mapping described in previous work for uncorrelated patterns, W = N−1/2ξ (Barra et al., 2012) (green); and (iv) the top p PCA components of the training data (pink). In Fig. 4(c), (d) we compare generated sample images from two RBMs, each with p = 50 hidden units but different initial weights (random in (c) and the HN mapping in (d)). The quality of samples in Fig. 4(d) reflect the efficient training of the Hopfield initialized RBM.
Fig. 4(a), (b) show that the Hopfield initialized weights provide an advantage over other approaches during the early stages of training. The PCA and Hebbian initializations start at much lower values of the objective and require one or more epochs of training to perform similarly (Fig. 4(a) inset), while the randomly initialized RBM takes > 25 epochs to reach a similar level. All initializations ultimately reach the same value. This is noteworthy because the proposed weight initialization is fast compared to conventional RBM training. PCA performs best for intermediate training times.
Despite being a common choice, the random initialization trains surprisingly slowly, taking roughly 40 epochs in Fig. 4(a), and in Fig. 4(b) we had to increase the basal learning rate η0 = 10−4 by a factor of 5 for the first 25 epochs due to slow training. The non-random initializations, by comparison, arrive at the same maximum value much sooner. The relatively small change over training for the Hopfield initialized weights supports the idea that they may be near a local optimum of the objective, and that conventional training may simply be mildly tuning them (Fig. 3).
That the HN initialization performs well at 0 epochs suggests that the p Hopfield patterns concisely summarize the dataset. This is intuitive, as the projection rule encodes the patterns (and nearby states) as high probability basins in the free energy landscape of Eq. (1). As the data itself is clustered near the patterns, these basins should model the true data distribution well. Overall, our results suggest that the HN correspondence provides a useful initialization for generative modelling with binary-gaussian RBMs, displaying excellent performance with minimal training.

heading: 3.2 CLASSIFICATION OBJECTIVE
text: As with the generative objective, we find that the Hopfield initialization provides an advantage for classification tasks. Here we consider the closely related MNIST classification problem. The goal is to train a model on the MNIST Training dataset which accurately predicts the class of presented images. The key statistic is the number of misclassified images on the MNIST Testing dataset.
We found relatively poor classification results with the single (large) RBM architecture from the preceding Section 3.1. Instead, we use a minimal product-of-experts (PoE) architecture as described in Hinton (2002): the input data is first passed to 10 RBMs, one for each class µ. This “layer of RBMs” functions as a pre-processing layer which maps the high dimensional sample s to a feature vector f(s) ∈ R10. This feature vector is then passed to a logistic regression layer in order to predict the class of s. The RBM layer and the classification layer are trained separately.
The first step is to train the RBM layer to produce useful features for classification. As in Hinton (2002), each small RBM is trained to model the distribution of samples from a specific digit class µ. We use CD-20 generative training as in Section 3.1, with the caveat that each expert is trained solely on examples from their respective class. Each RBM connects N binary visible units to k gaussian hidden units, and becomes an “expert” at generating samples from one class. To focus on the effect of interlayer weight initialization, we set the layer biases to 0.
After generative training, each expert should have relatively high probability p(µ)θ (sa) for sample digits sa of the corresponding class µ, and lower probability for digits from other classes. This idea is used to define 10 features, one from each expert, based on the log-probability of a given sample
under each expert, ln p(µ)θ (sa) = −βH(µ)(sa)− lnZ(µ). Note that β and lnZ(µ) are constants with respect to the data and thus irrelevant for classification. For a binary-gaussian RBM, H(µ)(sa) has the simple form Eq. (10), so the features we use are
f (µ)(s) = ∑ i ∑ j ∑ ν W (µ) iν W (µ) jν sisj = ||s TW (µ)||2. (16)
With the feature map defined, we then train a standard logistic regression classifier (using scikitlearn (Pedregosa et al., 2011)) on these features. In Fig. 5, we report the classification error on the MNIST Testing set of 10,000 images (held out during both generative and classification training). Note the size p = 10 of the feature vector is independent of the hidden dimension k of each RBM, so the classifier is very efficient.
Despite this relatively simple approach, the PoE initialized using the orthogonalized Hopfield patterns (“Hopfield PoE”) performs fairly well (Fig. 5, blue curves), especially as the number of subpatterns is increased. We found that generative training beyond 50 epochs did not significantly improve performance for the projection HN or PCA. (in Fig. E.1, we train to 100 epochs and also display the aforementioned “Hebbian” initial condition, which performs much worse for classification). Intuitively, increasing the number of hidden units k increases classification performance independent of weight initialization (with sufficient training).
For k fixed, the Hopfield initialization provides a significant benefit to classification performance compared to the randomly initialized weights (purple curves). For few sub-patterns (circles k = 10 and squares k = 20), the Hopfield initialized models perform best without additional training and until 1 epoch, after which PCA (pink) performs better. When each RBM has k = 100 hidden features (triangles), the Hopfield and PCA PoE reach 3.0% training error, whereas the randomly initialized PoE reaches 4.5%. However, the Hopfield PoE performs much better than PCA with minimal training, and maintains its advantage until 10 epochs, after which they are similar. Interestingly, both the Hopfield and PCA initialized PoE with just k = 10 encoded patterns performs better than or equal to the k = 100 randomly initialized PoE at each epoch despite having an order of magnitude fewer trainable parameters. Finally, without any generative training (0 epochs), the k = 100 Hopfield PoE performs slightly better (4.4%) than the k = 100 randomly initialized PoE with 50 epochs of training.

heading: 4 DISCUSSION
text: We have presented an explicit, exact mapping from projection rule Hopfield networks to Restricted Boltzmann Machines with binary visible units and gaussian hidden units. This provides a generalization of previous results which considered uncorrelated patterns (Barra et al., 2012), or special cases of correlated patterns (Agliari et al., 2013; Mézard, 2017). We provide an initial characterization of the reverse map from RBMs to HNs, along with a matrix-factorization approach to construct approximate associated HNs when the exact reverse map is not possible. Importantly, our HN to
RBM mapping can be applied to correlated patterns such as those found in real world datasets. As a result, we are able to conduct experiments (Section 3) on the MNIST dataset which suggest the mapping provides several advantages.
The conversion of an HN to an equivalent RBM has practical utility: it trades simplicity of presentation for faster processing. The weight matrix of the RBM is potentially much smaller than the HN (Np elements instead ofN(N −1)/2). More importantly, proper sampling of stochastic trajectories in HNs requires asynchronous updates of the units, whereas RBM dynamics can be simulated in a parallelizable, layer-wise fashion. We also utilized the mapping to efficiently estimate the partition function of the Hopfield network (Fig. 2) by summing out the spins after representing it as an RBM.
This mapping also has another practical utility. When used as an RBM weight initialization, the HN correspondence enables efficient training of generative models (Section 3.1, Fig. 4). RBMs initialized with random weights and trained for a moderate amount of time perform worse than RBMs initialized to orthogonalized Hopfield patterns and not trained at all. Further, with mild training of just a few epochs, Hopfield RBMs outperform conventionally initialized RBMs trained several times longer. The revealed initialization also shows advantages over alternative non-random initializations (PCA and the “Hebbian” Hopfield mapping) during early training. By leveraging this advantage for generative tasks, we show that the correspondence can also be used to improve classification performance (Section 3.2, Fig. 5, Appendix E.3).
Overall, the RBM initialization revealed by the mapping allows for smaller models which perform better despite shorter training time (for instance, using fewer hidden units to achieve similar classification performance). Reducing the size and training time of models is critical, as more realistic datasets (e.g. gene expression data from single-cell RNA sequencing) may require orders of magnitude more visible units. For generative modelling of such high dimensional data, our proposed weight initialization based on orthogonalized Hopfield patterns could be of practical use. Our theory and experiments are a proof-of-principle; if they can be extended to the large family of deep architectures which are built upon RBMs, such as deep belief networks (Hinton et al., 2006) and deep Boltzmann machines (Salakhutdinov & Hinton, 2009), it would be of great benefit. This will be explored in future work.
More broadly, exposing the relationship between RBMs and their representative HNs helps to address the infamous interpretability problem of machine learning which criticizes trained models as “black boxes”. HNs are relatively transparent models, where the role of the patterns as robust dynamical attractors is theoretically well-understood. We believe this correspondence, along with future work to further characterize the reverse map, will be especially fruitful for explaining the performance of deep architectures constructed from RBMs.

heading: ACKNOWLEDGMENTS
text: The authors thank Duncan Kirby and Jeremy Rothschild for helpful comments and discussions. This work is supported by the National Science and Engineering Research Council of Canada (NSERC) through Discovery Grant RGPIN 402591 to A.Z. and CGS-D Graduate Fellowship to M.S.

heading: A HOPFIELD NETWORK DETAILS
text: Consider p < N N -dimensional binary patterns {ξµ}pµ=1 that are to be “stored”. From them, construct the N ×p matrix ξ whose columns are the p patterns. If they are mutually orthogonal (e.g. randomly sampled patterns in the large N → ∞ limit), then choosing interactions according to the Hebbian rule, JHebb = 1N ξξ T , guarantees that they will all be stable minima of H(s) = − 12s TJs, provided α ≡ p/N < αc, where αc ≈ 0.14 (Amit et al., 1985). If they are not mutually orthogonal (referred to as correlated), then using the “projection rule” JProj = ξ(ξT ξ)−1ξT guarantees that they will all be stable minima of H(s), provided p < N (Kanter & Sompolinsky, 1987; Personnaz et al., 1986). Note JProj → JHebb in the limit of orthogonal patterns. In the main text, we use J as shorthand for JProj.
We provide some relevant notation from Kanter & Sompolinsky (1987). Define the overlap of a state s with the p patterns as m(s) ≡ 1N ξ
Ts, and define the projection of a state s onto the p patterns as a(s) ≡ (ξT ξ)−1ξTs ≡ A−1m(s). NoteA ≡ ξT ξ is the overlap matrix, andmµ,aµ ∈ [−1, 1]. We can re-write the projection rule Hamiltonian Eq. (1) as
H(s) = −N2m(s) · a(s). (A.1)
For simplicity, we include the self-interactions rather than keeping track of their omission; the results are the same as N → ∞. From Eq. (A.1), several quadratic forms can be written depending on which variables one wants to work with:
I. H(s) = −N 2 2 m T (ξT ξ)−1m
II. H(s) = − 12a T (ξT ξ)a
These are the starting points for the alternative Boltzmann Machines (i.e. not RBMs) presented in Appendix C.

heading: B ADDITIONAL HN TO RBM MAPPINGS
text: We used QR factorization in the main text to establish the HN to RBM mapping. However, one can use any decomposition which satisfies
JProj = UU T (B.1)
such thatU ∈ RN×p is orthogonal (orthogonal for tall matrices meansUTU = I). In that case, U becomes the RBM weights. We provide two simple alternatives below, and show they are all part of the same family of orthogonal decompositions.
“Square root” decomposition: Define the matrixK ≡ ξ(ξT ξ)−1/2. Note thatK is orthogonal, and that JProj =KKT .
Singular value decomposition: More generally, consider the SVD of the pattern matrix ξ:
ξ = UΣV T (B.2)
where U ∈ RN×p, V ∈ Rp×p store the left and right singular vectors (respectively) of ξ as orthogonal columns, and Σ ∈ Rp×p is diagonal and contains the singular values of ξ. Note in the limit of orthogonal patterns, we have Σ = √ NI . This decomposition gives several relations for quantities of interest: A ≡ ξT ξ = V Σ2V T
JHebb ≡ 1
N ξξT =
1
N UΣ2UT
JProj ≡ ξ(ξT ξ)−1ξT = UUT .
(B.3)
The last line is simply the diagonalization of JProj, and shows that our RBM mapping is preserved if we swap the Q from QR with U from SVD. However, since there are p degenerate eigenvalues σ2 = 1,U is not unique - any orthogonal basis for the 1-eigenspace can be chosen. ThusU ′ = UO where O is orthogonal is also valid, and the QR decomposition and “square root” decomposition correspond to particular choices ofO.

heading: C HN TO BM MAPPINGS USING ALTERNATIVE REPRESENTATIONS OF THE HOPFIELD HAMILTONIAN
text: In addition to the orthogonalized representation in the main text, there are two natural representations to consider based on the pattern overlaps and projections introduced in Appendix A. These lead to generalized Boltzmann Machines (BMs) consisting of N original binary spins and p continuous variables. These representations are not RBMs as the continuous variables interact with each other. We present them for completeness.
“Overlap” Boltzmann Machine: Writing H(s) = −N 2 2 m T (ξT ξ)−1m, we have
Z = ∑ {s} exp
( 1
2 (β √ Nm)T
( β
N ξT ξ
)−1 (β √ Nm) ) . (C.1)
Applying the gaussian integral identity, Z = √ det (ξT ξ) ∑ {s} ∫ e − β2N ∑ µ,ν(ξ T ξ)µνλµλν+ β√ N ∑ µ λµ ∑ i ξiµsi ∏ µ dλµ√ 2πN/β , (C.2)
from which we identify the BM Hamiltonian
H(s,λ) = 1
2N ∑ µ,ν (ξT ξ)µνλµλν − 1√ N ∑ µ ∑ i ξiµsiλµ. (C.3)
This is the analog of Eq. (6) in the main text for the “overlap” representation. Note we can also sum out the binary variables in Eq. (C.2), which allows for an analogous expression to Eq. (8),
F0({λµ}) = 1
2N ∑ µ,ν (ξT ξ)µνλµλν − 1 β ∑ i ln cosh ( β√ N ∑ µ ξiµλµ ) . (C.4)
Curiously, we note that we may perform a second Gaussian integral on Eq. (C.2), introducing new auxiliary variables τν to remove the interactions between the λµ variables:
Z = √ det ( β2
N ξT ξ )∑ {s} ∫ ∫ e − β2 τ T τ+ β√ N λT ξT s+i β√ N λT (ξT ξ)1/2τ ∏ ν dτν√ 2π ∏ µ dλµ√ 2π . (C.5)
Eq. (C.5) describes a three-layer RBM with complex interactions between the λµ and τν variables, a representation which could be useful in some contexts.
“Projection” Boltzmann Machine: Proceeding as above but for H(s) = − 12a T (ξT ξ)a, one finds
Z = det (ξT ξ) −1/2∑
{s}
∫ e− β 2 λ T (ξT ξ)−1λ+βλT (ξT ξ)−1ξT s ∏ µ dλµ√ 2π/β , (C.6)
which corresponds to the BM Hamiltonian
H(s,λ) = 1
2 λT (ξT ξ)−1λ− λT (ξT ξ)−1ξTs. (C.7)
The analogous expression to Eq. (8) in this case is
F0(λ) = 1 2 λT (ξT ξ)−1λ− 1 β ∑ i ln cosh ( β[ξ(ξT ξ)−1λ]i ) . (C.8)

heading: D RBM TO HN DETAILS
text: D.1 INTEGRATING OUT THE HIDDEN VARIABLES
The explanation from Mehta et al. (2019) for integrating out the hidden variables of an RBM is presented here for completeness. For a given binary-gaussian RBM defined by HRBM(s,λ) (as in Eq. (9)), we have p(s) = Z−1 ∫ e−βHRBM(s,λ)dλ. Consider also that p(s) = Z−1e−βH(s) for some unknown H(s). Equating these expressions gives,
H(s) = − ∑ i bisi − 1 β ∑ µ ln (∫ e−β 1 2λ 2 µ+β ∑ iWiµsiλµdλµ ) . (D.1)
Decompose the argument of ln(·) in Eq. (D.1) by defining qµ(λµ), a gaussian with zero mean and variance β−1. Writing tµ = β ∑ iWiµsi, one observes that the second term (up to a constant) is a sum of cumulant generating functions, i.e.
Kµ(tµ) ≡ ln 〈etµλµ〉qµ = ln (∫ qµe tµλµdλµ ) . (D.2)
These can be written as a cumulant expansion, Kµ(tµ) = ∑ n=1 κ (n) µ tnµ n! , where κ (n) µ = ∂Kµ ∂tµ |tµ=0 is the nth cumulant of qµ. However, since qµ(λµ) is gaussian, only the second term remains, leaving
Kµ(tµ) = 1 β
(β ∑ iWiµsi) 2
2 . Putting this all together, we have
H(s) = − ∑ i bisi − 1 2 ∑ i ∑ j ∑ µ WiµWjµsisj , (D.3)
Note that in general, qµ(λµ) need not be gaussian, in which case the resultant Hamiltonian H(s) can have higher order interactions (expressed via the cumulant expansion).
D.2 APPROXIMATE REVERSE MAPPING
Suppose one has a trial solutionBp =WX to the approximate binarization problem Eq. (12), with error matrix E ≡ Bp − sgn(Bp). We consider two cases depending on ifW is orthogonal.
Case 1: If W is orthogonal, then starting from Eq. (11) and applying Eq. (13), we have J = WW T = Bp(X TX)−1BTp . Using I =W TW , we get
J = Bp(B T p Bp) −1BTp . (D.4)
Thus, the interactions between the visible units is the familiar projection rule used to store the approximately binary patternsBp. “Storage” means the patterns are stable fixed points of the deterministic update rule st+1 ≡ sgn(Jst). We cannot initialize the network to a non-binary state. Therefore, the columns ofB = sgn(Bp) are our candidate patterns. To test if they are fixed points, consider
sgn(JB) = sgn(JBp − JE) = sgn(Bp − JE).
We need the errorE to be such that JE will not alter the sign ofBp. Two sufficient conditions are:
(a) small error: |(JE)iµ| < |(Bp)iµ|, or (b) error with compatible sign: (JE)iµ(Bp)iµ < 0.
When either of these conditions hold for each element, we have sgn(JB) = sgn(Bp) = B, so that the candidate patternsB are fixed points. It remains to show that they are also stable (i.e. minima).
Case 2: If W is not orthogonal but its singular values remain close to one, then Löwdin Orthogonalization (also known as Symmetric Orthogonalization) (Löwdin, 1970) provides a way to preserve the HN mapping from Case 1 above.
Consider the SVD of W : W = UΣV T . The closest matrix to W (in terms of Frobenius norm) with orthogonal columns is L = UV T , and the approximation W ≈ L is called the Löwdin Orthogonalization of W . Note the approximation becomes exact when all the singular values are one. We then writeWW T ≈ UUT , and the orthogonalW approach of Case 1 can then be applied. On the other hand, W may be strongly not orthogonal (singular values far from one). If it is still full rank, then its pseudo-inverse W † = (W TW )−1W T = V Σ−1UT is well-defined. Repeating the steps from the orthogonal case, we note here XTX = BTp (W
†)TW †Bp. Defining C ≡ (W †)TW † = UΣ−2UT , we arrive at the corresponding result,
J = Bp(B T p CBp) −1BTp . (D.5)
This is analogous to the projection rule but with a “correction factor” C. However, it is not immediately clear how C affects pattern storage. Given the resemblance between JProj and Eq. (D.4) (relative to Eq. (D.5)), we expect that RBMs trained with an orthogonality constraint on the weights may be more readily mapped to HNs.
D.3 EXAMPLE OF THE APPROXIMATE REVERSE MAPPING
In the main text we introduced the approximate binarization problem Eq. (12), the solutions of which provide approximately binary patterns through Eq. (13). To numerically solve Eq. (12) and obtain a candidate solutionX∗, we perform gradient descent on a differentiable variant.
Specifically, we replace sgn(u) with tanh(αu) for large α. Define E = WX − tanh(αWX) as in the main text. Then the derivative of the “softened” Eq. (12) with respect toX is
G(X) = 2W T (E − αE sech2(WX))). (D.6)
Given an initial conditionX0, we apply the update rule
Xt+1 =Xt − γG(Xt) (D.7)
until convergence to a local minimumX∗.
In the absence of prior information, we consider randomly initialized X0. Our preliminary experiments using Eq. (D.7) to binarize arbitrary RBM weightsW have generally led to high binarization error. This is due in part to the difficulty in choosing a good initial condition X0, which will be explored in future work.
To avoid this issue, we consider the case of a Hopfield-initialized RBM following CD-k training. At epoch 0, we in fact have the exact binarization solution X∗ = R (from the QR decomposition, Eq. (2)), which recovers the encoded binary patterns ξ. We may use X0 = R, as an informed initial condition to Eq. (D.7), to approximately binarize the weight at later epochs and monitor how the learned patterns change.
In Fig. D.1, we give an example of this approximate reverse mapping for a Hopfield-initialized RBM following generative training (Fig. 4). Fig. D.1(a) shows the p = 10 encoded binary patterns, denoted below by ξ0, and Fig. D.1(b) shows the approximate reverse mapping applied to the RBM weights at epoch 10. We denote these nearly binary patterns by ξ10. Interestingly, some of the non-binary regions in Fig. D.1(b) coincide with features that distinguish the respective pattern. For example, the strongly “off” area in the top-right of the “six” pattern.
We also considered the associative memory performance of the projection HNs built from the the patterns ξ0, ξ10 from Fig. D.1. Specifically, we test the extent to which images from the MNIST Test dataset are attracted to the correct patterns. Ideally, each pattern should attract all test images of the corresponding class. The results are displayed in Fig. D.2 and elaborated in the caption.
The results in Fig. D.2 suggest that p = 10 patterns may be too crude for the associative memory network to be used as an accurate MNIST classifier (as compared to e.g. Fig. 5). Notably, the HN constructed from ξ10 performs about 3% worse than the HN constructed from ξ0, although this performance might improve with a more sophisticated method for the associative memory task. There may therefore be a cost, in terms of associative memory performance, to increasing the generative functionality of such models (Fig. 4). Our results from Appendix D.2 indicate that incorporating an orthogonality constraint on the weights during CD-k generative training may provide a way to preserve or increase the associative memory functionality. This will be explored in future work.

heading: E RBM TRAINING
text: Consider a general binary-gaussian RBM with N visible and p hidden units. The energy function is
H(s,λ) = 1
2 ∑ µ (λµ − cµ)2 − ∑ i bisi − ∑ µ ∑ i Wiµsiλµ. (E.1)
First, we note the Gibbs-block update distributions for sampling one layer of a binary-gaussian RBM given the other (see e.g. Melchior et al. (2017)),
Visible units: p(si = 1|λ) = 11+e−2βxi , where xi ≡ ∑ µWiµλµ + bi defines input to si,
Hidden units: p(λµ = λ|s) ∼ N (hµ, β−1), where hµ ≡ ∑ iWiµsi + cµ defines the input to λµ.
E.1 GENERATIVE TRAINING
For completeness, we re-derive the binary-gaussian RBM weight updates, along the lines of Melchior et al. (2017). We want to maximize L ≡ 1M ∑ a ln pθ(sa). The contribution for a single
datapoint sa has the form ln pθ(sa) = ln(C−1 ∫ e−βH(sa,λ)dλ)− lnZ with C ≡ (2π/β)p/2. The gradient with respect to the model is
∂ ln pθ(sa)
∂θ =
∫ (−β ∂H(sa,λ)∂θ )e −βH(sa,λ) ∏ µ dλµ∫
e−βH(sa,λ) ∏ µ dλµ
− β ∑ {s} ∫ siλµe −βH(sa,λ) ∏ µ dλµ
Z (E.2)
We are focused on the interlayer weights, with ∂H(s,λ)∂Wiµ = −siλµ, so
∂ ln pθ(sa)
∂θ = β(sa)i
∫ (λµe −βH(sa,λ) ∏ µ dλµ∫
e−βH(sa,λ) ∏ µ dλµ
− β ∑ {s} ∫ siλµe −βH(sa,λ) ∏ µ dλµ
Z = β(sa)i〈λµ|s = sa〉model − β〈siλµ〉model. (E.3)
The first term is straightforward to compute: 〈λµ|s = sa〉model = ∑ iWiµ(sa)i + cµ. The second term is intractable and needs to be approximated. We use contrastive divergence (Carreira-Perpinan & Hinton, 2005; Hinton, 2012): 〈siλµ〉model ≈ s(k)i λ (k) µ . Here k denotes CD-k – k steps of Gibbsblock updates (introduced above) – from which s(k)i , λ (k) µ comprise the final state. We evaluate both terms over mini-batches of the training data to arrive at the weight update rule Eq. (15).
E.2 GENERATIVE PERFORMANCE We are interested in estimating the objective function L ≡ 1M ∑ a ln pθ(sa) during training. As above, we split L into two terms,
L = ln ( C−1 ∫ e−βH(sa,λ)
∏ µ dλµ
) − lnZ, (E.4)
with C ≡ (2π/β)p/2. The first term evaluates to
β
M M∑ a=1 (bTsa + 1 2 ||c+W Tsa||2)− β 2 ||c||2, (E.5)
which can be computed deterministically. lnZ, on the other hand, needs to be estimated. For this we perform Annealed Importance Sampling (AIS) (Neal, 2001) on its continuous representation
Z = C−12N ∫ e− β 2 ∑ µ(λµ−cµ) 2+ ∑ i ln cosh(β( ∑ µWiµλµ+bi)) ∏ µ dλµ. (E.6)
For AIS we need to specify the target distribution’s un-normalized log-probabilities
ln(Zp(λ)) = N ln 2− p 2 ln
( 2π
β
) − β
2 ∑ µ (λµ − cµ)2 + ∑ i ln cosh
( β (∑ µ Wiµλµ + bi )) ,
(E.7) as well as an initial proposal distribution to anneal from, which we fix as a p-dimensional unit normal distribution N (0, I).
E.3 CLASSIFICATION PERFORMANCE
Here we provide extended data (Fig. E.1) on the classification performance shown in Fig. 5. As in the main text, color denotes initial condition (introduced in Section 3.1) and shape denotes the number of sub-patterns. Here we train for each 100 epochs, which allows convergence of most of the curves. We also include the “Hebbian” initialization (green curves).
Figure E.1: Product-of-experts classification performance (extended data).
Notably, the Hebbian initialization performs quite poorly in the classification task (as compared to direct generative objective, Fig. 4). In particular, for the 100 sub-patterns case, where the projection
HN performs best, the Hebbian curve trains very slowly (still not converged after 100 epochs) and lags behind even the 10 sub-pattern Hebbian curve for most of the training. This emphasizes the benefits of the projection rule HN over the Hebbian HN when the data is composed of correlated patterns, which applies to most real-world datasets.

Review summaries:
This paper shows a relationship between the project rule weights of a Hopfield network (HN) and the interaction weights in a corresponding restricted Boltzmann machine (RBM). The mapping from HN to RBM is facilitated by realising that the partition function of BN can be seen as the partition function of a binary-continuous (Bernoulli-Gaussian) RBM. The authors comments on the mapping from RBM to BN. The experiments show the advantages of training RBM with weights initialised from BN projection weights in generation and classification.