Random active path model of deep neural networks with diluted binary synapses

TL;DR

This paper introduces a random active path model for deep neural networks with binary synapses, analyzing how connection removal affects collective properties and phase transitions, providing insights into generalization performance.

Contribution

It proposes a novel random active path model to study deep neural networks with binary synapses under perturbation, revealing phase transition behavior.

Findings

Identifies a critical perturbation value separating spin glass and paramagnetic regimes.

Shows the paramagnetic phase likely has poor generalization.

Demonstrates a first-order phase transition in the model.

Abstract

Deep learning has become a powerful and popular tool for a variety of machine learning tasks. However, it is challenging to understand the mechanism of deep learning from a theoretical perspective. In this work, we propose a random active path model to study collective properties of deep neural networks with binary synapses, under the removal perturbation of connections between layers. In the model, the path from input to output is randomly activated, and the corresponding input unit constrains the weights along the path into the form of a $p$-weight interaction glass model. A critical value of the perturbation is observed to separate a spin glass regime from a paramagnetic regime, with the transition being of the first order. The paramagnetic phase is conjectured to have a poor generalization performance.

Equations

$$\mathbf{y}=f_{L-1}\left({\mathbf{W}}^{L-1}{f}_{L-2}({\mathbf{W}}^{L-2}\cdots f_1({\mathbf{W}}^1\mathbf{v}))\right),$$

$$y=\sum _{a=1}^{\Psi }{v}^a\prod _{k=1}^{L-1}{W}_a^k,$$

$$\mathcal{H}(\mathbf{W})=-\sum _{a=1}^{n^3}{A}_a{v}^a\prod _{i\in \partial a}{W}_i,$$

$$P(E_{1},E_{2},\ldots,E_{\mathcal{N}})=\int_{-\infty}^{+\infty}\prod_{a=1}^{\mathcal{N}}\frac{d\hat{E}_{a}}{2\pi}\exp\left(\sum_{a}{\rm i}\hat{E}_{a}E_{a}+n^{p}\prod_{l}\xi_{l}(1-p_{l})\Bigl{(}\prod_{a}\cosh({\rm i}\hat{E}_{a})-1\Bigr{)}\right),$$

$$\begin{split}Z^{new}&=\sum_{W_{i}}\sum_{\mathbf{W}}\exp\bigl{(}{\sum_{a=1}^{M}\beta v^{a}\prod_{k\in\partial a}W_{k}+\beta\sum_{b=1}^{C}v^{b}W_{i}\prod_{j\in\partial b\backslash i}W_{j}}\bigr{)}\\ &=Z^{old}\sum_{W_{i}}\prod_{b}\sum_{\{W_{j}\}:j\in\partial b\backslash i}\prod_{j\in\partial b\backslash i}\biggl{[}\frac{1+W_{j}m_{j\rightarrow b}}{2}\biggr{]}\cdot e^{\beta v^{b}W_{i}\prod_{j\in\partial b\backslash i}W_{j}}\\ &=Z^{old}\biggl{[}\prod_{b}\bigl{[}\cosh\beta v^{b}\bigl{(}1+\tanh\beta v^{b}\prod_{j\in\partial b\backslash i}m_{j\rightarrow b}\bigr{)}\bigr{]}+\prod_{b}\bigl{[}\cosh\beta v^{b}\bigl{(}1-\tanh\beta v^{b}\prod_{j\in\partial b\backslash i}m_{j\rightarrow b}\bigr{)}\bigr{]}\biggr{]},\end{split}$$

$$-\beta\Delta F_{i}=\ln\frac{Z^{new}}{Z^{old}}=\ln\biggl{[}\prod_{b\in\partial i}\bigl{[}\cosh\beta v^{b}(1+\hat{m}_{b\rightarrow i})\bigr{]}+\prod_{b\in\partial i}\bigl{[}\cosh\beta v^{b}(1-\hat{m}_{b\rightarrow i})\bigr{]}\biggr{]}.$$

$$\begin{split}Z^{new}&=\sum_{\mathbf{W}}\exp\bigl{(}{\beta\sum_{b=1}^{M}v^{b}\prod_{k\in\partial b}W_{k}+\beta v^{a}\prod_{i\in\partial a}W_{i}}\bigr{)}\\ &=\sum_{\mathbf{W}}e^{\beta\sum_{b=1}^{M}v^{b}\prod_{k\in\partial b}W_{k}}\sum_{\mathbf{W}}\frac{e^{\beta\sum_{b=1}^{M}v^{b}\prod_{k\in\partial b}W_{k}}}{\sum_{\mathbf{W}}e^{\beta\sum_{b=1}^{M}v^{b}\prod_{k\in\partial b}W_{k}}}e^{\beta v^{a}\prod_{i\in\partial a}W_{i}}\\ &=Z^{old}\sum_{\mathbf{W}}P(\mathbf{W})e^{\beta v^{a}\prod_{i\in\partial a}W_{i}}\\ &=Z^{old}\sum_{\{W_{i}\}:i\in\partial a}\prod_{i\in\partial a}\biggl{[}\frac{1+m_{i\rightarrow a}W_{i}}{2}\biggr{]}e^{\beta v^{a}\prod_{i\in\partial a}W_{i}}\\ &=Z^{old}\cdot\cosh\beta v^{a}\bigl{(}1+\tanh\beta v^{a}\prod_{i\in\partial a}m_{i\rightarrow a}\bigr{)}.\end{split}$$

$$F\equiv -\ln Z=\sum _i\Delta F_i-\sum _a(|\partial a|-1)\Delta F_a.$$

$$m_{i\to a}$$

$${\hat{m}}_{b\to i}$$

$$E$$

$$\Delta E_i$$

$$\Delta E_a$$

$$\displaystyle\begin{split}\mathcal{G}_{i}(x)&=\sum_{b\in\partial i}\Biggl{[}v^{b}\sinh\beta v^{b}(1+x\hat{m}_{b\rightarrow i})+xv^{b}\cosh\beta v^{b}(1-\tanh^{2}\beta v^{b})\\ &\times\prod_{j\in\partial b\backslash i}m_{j\rightarrow b}\Biggr{]}\prod_{a\in\partial i\backslash b}\cosh\beta v^{a}(1+x\hat{m}_{a\rightarrow i}),\end{split}$$

$${\mathcal{H}}_i(x)$$

Full text

I Introduction

Deep learning (e.g., with many layers of neural networks) works very well in areas from speech recognition, image classification, to drug discovery, medical image analysis, particle discovery, automatic game playing, and many others Goodfellow et al. (2016). This is due to the available large dataset for training and efficient hardware design such as GPU to accelerate training, rather than breakthrough in theoretical foundations. Theoretical efforts are recently devoted to address partially the origin of these impressive performances of deep networks, relying on simple assumptions Choromanska et al. (2014); Chaudhari and Soatto (2015); Kawaguchi (2016); Soudry and Carmon (2016); Vardan et al. (2016); Poggio et al. (2018); Poole et al. (2016).

One of promising perspectives looking at deep learning is the concept of redundancy, i.e., the deep architecture is robust, to some extent, in the absence of a fraction of connections between layers. In other words, the generalization ability measured by the output error on test (unseen) dataset does not significantly change until a sufficient amount of connections are removed. Indeed, it was recently revealed that deep convolutional neural networks mimicking the ventral visual pathway are robust to a number of weight perturbations in the higher convolutional layers Cheney et al. (2017). We also observed that a multi-layer perceptron network classifying handwritten digits has such redundancy behavior (Supplemental Material). This suggests that the redundancy property may be a general principle of deep computation. In terms of synaptic activities, this kind of redundancy is expected to appear from interactions between synapses, in a supervised learning system.

From a theoretical perspective, an interesting question arises, i.e., how the statistical property of a deep neural network changes with respect to the removal perturbation of synapses between adjacent layers. We address this question by proposing a simple random active path (RAP) model, in which we construct randomly and independently each active path from the input at the bottom layer to the output at the top layer ($L$ layers in total), and weights along each path are constrained by the corresponding random input of that path. The RAP model is thus formulated as a $p$-weight interaction model ($p=L-1$ refers to the depth of the deep network). For simplicity, we assume that the path is randomly activated in the sense that units in each layer are randomly activated with a layer-dependent activation probability and the activation is thus not input-dependent. The activation is defined in terms of ReLU function commonly used in deep learning Nair and Hinton (2010). In addition, the activation of a path also depends on the removal perturbation of synapses that can make a path inactivated. One significant difference from previous $p$-spin glass models Gardner (1985); Kirkpatrick and Thirumalai (1987); Franz et al. (2001); Montanari, A. and Ricci-Tersenghi, F. (2003) is that the degree (the number of paths one weight is involved in) distribution of weights in the model depends on both the layer-dependent activation probability and the perturbation level of deep neural networks, such that there are multiple peaks in the distribution, and the degree can vary in a wide range of values.

The RAP model can be analytically studied by mean-field methods, revealing a critical value of the perturbation separating a paramagnetic phase from a spin glass phase. The paramagnetic phase is conjectured to be related to the non-robust regime of deep neural networks under weight perturbations. This result may also explain the performance of the dropconnect algorithm Wan et al. (2013), where a finite fraction of connections is stochastically turned off during training, when different fractions are used.

The paper is structured as follows. We first introduce the RAP model to study collective properties of deep neural networks with binary synapses, by varying the perturbation level. Theoretical predictions match very well Monte-Carlo simulation results. We give concluding remarks and future perspectives in the final section.

II A random active path model

In this section, we propose a random active path model to qualitatively understand the collective behavior of deep neural networks with diluted binary synapses, i.e., to address how the statistical property of a deep neural network changes with respect to the removal perturbation of synapses between adjacent layers. We consider a deep network model with $L=4$ layers of fully-connected feedforward architecture. Each layer has $n_{l}(l=0,\ldots,L-1)$ units (so-called width of that layer). The deep network architecture is specified by $n_{0}$-$n_{1}$-$n_{2}$-$n_{3}$. The input is defined as an $n_{0}$-dimensional vector $\mathbf{v}$, and the weight matrix $\mathbf{W}^{l}$ with binary ($\pm 1$) entries specifies connections between layer $l$ and layer $l-1$. A bias can be incorporated into the weight matrix by assuming an additional constant input unit. The output at the final layer is expressed as:

[TABLE]

where $f_{l}(\cdot)$ is an element-wise ReLU function for units at layer $l$, defined as $f_{l}(u_{i})={\rm max}(0,u_{i})$ where $u_{i}$ is the weighted-sum input to the unit $i$ at layer $l$. For simplicity, we assume that all layers except the top one use ReLU transfer function, and the top layer has only one output unit that uses an identity map Kawaguchi (2016).

An active path refers to the path from one input unit to one output unit with the property that each connection on the path is present (each synapse is deleted with a dilution probability $p_{l}$, where $l$ specifies the weight population between two consecutive layers ($l$ and $(l-1)$-layer)) and all units along the path are activated because of ReLU activation. It follows that the network output defined in Eq. (1) can be re-expressed as Choromanska et al. (2014); Chaudhari and Soatto (2015)

[TABLE]

where $\Psi$ denotes the total number of active paths converged to the output unit, $v^{a}$ denotes the input specified to the $a$-th path for the output, and $W^{k}_{a}\in\{\pm 1\}$ denotes the entry of $\mathbf{W}^{k}$ that is used in the $a$-th path for the output.

A random construction of active paths is shown in Fig. 1, where only active paths are used to construct the Hamiltonian of the model composed of many constraints (represented by square nodes), and each of these constraints (e.g., $a$) is given by the input $v^{a}$ specified to the corresponding active path. In this construction process, we assume that activation of paths is independent of input Choromanska et al. (2014), realized by the units of each hidden layer activated independently with an activation probability, say $\xi_{l}$ that is layer-dependent and can be empirically estimated from deep neural network trainings (Fig. 2). The Hamiltonian of the RAP model can thus be written as

[TABLE]

where $n$ denotes the layer width (the same for all layers except the last one), $A_{a}$ is a binary value indicating the activation ($A_{a}=1$) or silent ($A_{a}=0$), and $\partial a$ denotes the set of the weights involved in the active path $a$. According to the random construction, the probability of a constructed active path is defined by $P(A_{a}=1)=\prod_{l}\xi_{l}(1-p_{l})$, where $\xi_{l}$ and $p_{l}$ are aforementioned activation and dilution probabilities, respectively. Note that the weight configuration $\mathbf{W}$ in the model is a subset of $\{\mathbf{W}^{l}\}$ in Eq. (1) unless all $n^{3}$ paths of the deep network become active. Here, we assume that $v^{a}$ are randomly selected from components of a random input vector $\mathbf{v}$ whose components follow independently a binary distribution ($P(v_{i}=\pm 1)=0.5$). During the random construction, some constraints may share the same value of $v^{a}$, i.e., the corresponding paths share a common input unit.

To build a naive relationship between the Hamiltonian we studied here and the loss function (e.g., absolute loss) used in training (Fig. 2), we assume that the true label is denoted by $Y_{t}=\pm\Lambda$ where $\Lambda>0$ is a predefined maximal output. We then define $\mathcal{H}=-{\rm sgn}(Y_{t})y$. It is easy to verify that the learning process minimizing the absolute loss ($\mathcal{C}=|Y_{t}-y|$) between the target and the actual output $y$ is equivalent to finding the minimal value of $\mathcal{H}$ in the RAP model. If we assume $Y_{t}$ is random, ${\rm sgn}(Y_{t})$ can be absorbed into the input and the model is statistically invariant. The random hinge loss can also be similarly analyzed.

Although it is challenging to prove the assumptions we made above reasonable in practical deep networks, the above $p$-weight interaction model still provides us a nice starting point to qualitatively understand complicated properties of deep neural networks. We will show interesting (non-trivial) properties of this simple model below.

We first evaluate the empirical values of the layer-dependent activation probability ($\xi_{l}$) from training a deep neural network with the same structure as our model (Fig. 1). The only difference is that, the binary target at the top layer is generated by a teacher network with the same structure fed with the corresponding input. The result is shown in Fig. 2, which implies that the first hidden layer has a higher activation probability (nearly independent of trials) than the subsequent layer. The activation probability in the subsequent layer varies in a relatively large range, suggesting a sparse distributed representation for efficient computation in the deep ReLU network Glorot et al. (2011). Therefore, in our model analysis, we fix $\xi_{1}=0.5$, and $\xi_{2}=0.05,0.1,0.15$. In fact, the value of $\xi_{1}=0.5$ can be derived by assuming a normal distribution for the weighted-sum $u_{i}$. The model defined in Eq. (3) can be expressed as a factor graph like the one shown in Fig. 1. The degree of a weight is thereby defined as the number of constraints the weight is involved in. The degree distribution of weights in a typical example is shown in Fig. 3 (a), which suggests that there exists three peaks corresponding to weights drawn from three different weight populations. The well-structured degree distribution is a feature of the RAP model constructed to mimic the behavior of a deep neural network.

The RAP model (Eq. (3)) is thus different from previous $p$-spin glass models Gardner (1985); Kirkpatrick and Thirumalai (1987); Franz et al. (2001); Montanari, A. and Ricci-Tersenghi, F. (2003) in three aspects. First, the model is interpreted in terms of active-path decomposition of a deep network’s output that is directly connected to training/test errors. In particular, the binary activation indicator $A_{a}$ is controlled by two deep-network-dependent parameters—the dilution probability and layer-dependent activation probability of units. Second, when an active path is randomly constructed, the quenched disorder comes from the corresponding input unit, thereby being physically explained. Lastly, the constructed model has a distinct three-peak degree distribution in the corresponding graphical representation, and each peak reflects the feature of a different layer.

The mean-field model defined in Eq. (3) can be solved by the cavity method Mézard and Parisi (2001); Huang (2017). The technical details are summarized in Appendix A. We consider the case of varying $p_{1}$ with zero $p_{2}$ and $p_{3}$. Because of layer-dependent activation probability, the quantitative behavior depends on which weight population is diluted. As shown in Fig. 3 (b) and (c), given the inverse temperature $\beta=1$, the RAP model has a paramagnetic phase when $p_{1}$ is larger than $0.9$, and the energy keeps decreasing as $p_{1}$ decreases. The theoretical prediction of the cavity method is in agreement with the simple Monte-Carlo simulations of the model (simulation details are given in Appendix B), despite deviations around the transition point. The deviations may be caused by the nature of the transition and sampling accuracy of the quenched Monte-Carlo procedure. First, the transition is of the first-order (Fig. 3 (c)), because there exists metastability of states (paramagnetic state versus spin glass state in which weights prefer to polarize towards a particular direction, but a spatial average over all magnetizations vanishes). Therefore, even in paramagnetic phase, depending on the details of the evolution, the quenched Monte-Carlo dynamics has a certain probability to be attracted by the spin glass state, which seems to have a lower energy than the paramagnetic state yet a weak contribution to the equilibrium measure in the paramagnetic phase (in terms of free energy taking into account the entropy contribution). Second, the accuracy of the Monte-Carlo sampling is affected by the number of collected configurations and relaxation steps, and how abruptly the energy landscape is changed. The effect is expected to be stronger around the transition point than the regime where one phase already overwhelms the other, which is reflected by the large error bars in Fig. 3 (c).

As shown in the inset of Fig. 3 (b), the entropy increases when more and more weights are pruned. In this case of large $p_{1}$, the deep network system becomes less constrained with a large number of candidate weight configurations, but the energy is relatively high. The entropy curve also displays a slight jump (about $0.03$ when $p_{1}$ increases from $0.9$ to $0.91$), which is a characteristic of the first-order phase transition. Overall, the equilibrium properties of the RAP model demonstrate that reducing the dilution probability will trigger a phase transition towards a spin glass state, and the paramagnetic phase located at a high dilution probability is clearly not favored, because it does not support the generalization ability of the network (one does not expect a randomly polarized (with equal probability) yet equilibrium configuration of weights yields a good generalization), and thus may connect to the high test errors when a sufficiently large fraction of connections are removed (Supplemental Material).

Finally, we explore the effects of different inverse-temperatures and intermediate-layer activation probabilities on the phase transition in Fig. 4. When the activation probability decreases, the critical dilution probability decreases as well. This implies that the regime of the paramagnetic phase expands, impairing the generalization performance of the network when an intermediate number of weights are deleted. On the other hand, the inverse-temperature is a measure of the stochastic noise level in learning process. A high temperature also expands the paramagnetic phase, thereby impairing the network’s robustness.

We further investigate the effects of inverse-temperatures on the validity of our results. With decreasing the temperature, the iteration of mean-field equations (Appendix A) will lose its stability (Fig. 5 (a)), for example, exhibiting an oscillatory behavior. In this case, the weight space may develop complex structures, i.e., the low-energy configurations form exponentially many clusters, and they may further condense onto a finite number of states, especially when the temperature is further lowered down. Then we should adopt the replica symmetry breaking assumption to study this complex situation (Appendix A). However, within the temperature range we explore in Fig. 3 and Fig. 4, the iteration is stable, leading to a converged value for the order parameter. Moreover, the entropy does not become negative (Fig. 5 (b)). The energy level predicted by the mean-field method coincides well with that reached by the Monte-Carlo dynamics, although at a slightly low temperature, the mean-field method, compared with the Monte-Carlo dynamics, can reach a lower energy level (the inset of Fig. 5 (b)). We also perform simulated annealing experiments on the constructed model to study the effects of cooling rates on the final reached energy levels (technical details are given in Appendix B). Interestingly, when the cooling becomes slower, the energy level the dynamics reaches does not show a significant change (Fig. 5 (a)), indicating that a dynamical ergodicity-breaking transition is not apparent.

In this work, we only consider $4$-layer feedforward neural networks, which correspond to a three-body interaction spin model. It is interesting to study the infinite depth limit, because in a standard $p$-spin interaction model where the coupling between spins follows a Gaussian distribution, and all possible coupling among $p$ spins are considered, Derrida showed that in the large-$p$ limit, the model behaves like a random energy model, in which the joint energy level distribution becomes factorized, and a frozen phase can be identified when the temperature is lowered down Derrida (1981). In our current context, according to the techniques used in Refs. Derrida (1981); Dominicis and Mottishaw (1987), the joint distribution of $\mathcal{N}=2^{N_{w}}$ energy levels can be expressed as

[TABLE]

where $n^{p}\prod_{l}\xi_{l}(1-p_{l})$ is exactly the number of active paths in the model, and the number is expected to be finite when $p$ tends to be infinite. It follows that the single energy level follows a Gaussian distribution with a fluctuation of the order of $[n^{p}\prod_{l}\xi_{l}(1-p_{l})]^{1/2}$ around zero, to leading order. To derive Eq. (4), we assume that the overlap $q$ between two configurations is smaller than one in the magnitude, therefore $q^{p}$ is a negligible term. We also use the gauge transformation Gross and Mezard (1984). It is not immediate to conclude that the energy levels are independent random variables from Eq. (4). Studying the neural network in the infinite depth limit is very interesting. In practice, training a very deep network is not easy due to the gradient vanishing problem. However, in the infinite depth limit, the energy levels may be organized into a non-trivial structure (e.g., a hierarchical organization like the one in glassy systems). The detailed exploration of this limit is beyond the current scope, and we leave that for future works.

III Discussion

The deep architecture is expected to be robust in the absence of a small fraction of connections between layers. In practice, not all connections (synapses) are used for transforming the original input, given that the number of parameters in modern deep learning architectures is typically much larger than the amount of training data. Given the redundancy property, the deep network has the chance of repairing and modification of subnetworks, without affecting the overall function. It was recently revealed that deep convolutional neural networks mimicking the ventral visual pathway are robust to a number of weight perturbations in the higher convolutional layers Cheney et al. (2017). This suggests that the redundancy property may be a general principle of deep computation.

Motivated by these empirical observations, we propose the RAP model to understand how the collective properties of deep neural networks change with the removal perturbation. Reducing the number of active paths drives the network to a paramagnetic phase where the weights behave like a noise with zero mean. In the opposite direction, the network enters a spin glass regime through a first-order phase transition. In the spin glass phase, the weights polarize to their own particular directions, rather than fluctuating with zero mean. This regime is meaningful in the sense that the network learns the feature of the input-output relationship hidden in the training data, with a configuration of particularly-polarized weights, which may connect to the robustness property of deep networks. In contrast, the paramagnetic phase does not have such properties, since weights behave like a random guess (taking both directions with equal probability). In other words, the redundancy is completely removed. In addition, the nature of the transition implies that around the transition, the paramagnetic phase is still stable and competes with the spin glass phase, which predicts that the learning should be difficult when the hyper-parameter such as dilution probability is chosen to be around the transition point.

Our model may connect to the dropconnect technique Wan et al. (2013), one popular regularization technique used during training deep networks. The dropconnect was proposed to stochastically turn off a finite fraction of connections during training, providing a sampling of model ensemble in the course of the training, thereby reducing over-fitting effects. When the dropconnect probability (the fraction of connections kept during training) is very small, the performance gets worse rapidly as the probability further decreases Wan et al. (2013), which is consistent with the paramagnetic phase as explained above. This suggests that one should choose a relatively large dropconnect probability to implement a sampling of good model ensemble in terms of its robustness in the learning performance.

The RAP model still needs to be improved by considering how activation of paths relies on the (structured) input and its specific label, together with conditional dependence of units’ activity at one layer on their previous layer. Thus this work encourages further refined models that can have quantitative predictions of the collective behavior observed in practical applications of deep learning. On the other hand, when the learning noise is weak, which corresponds to using a low temperature to explore the model’s properties, the model may exhibit rich equilibrium/dynamical properties, for which replica symmetry breaking calculations or glassy dynamics analysis may provide theoretical insights. These insights may also clarify how the robustness property of a deep network under training is affected. We leave this for future works.

Appendix A Mean field method to solve the RAP model

We present self-consistent mean-field equations in this appendix to analyze the statistical properties of the RAP model (Eq. (3)). These equations can be derived using the standard cavity method Mézard and Montanari (2009); Huang (2017). First, the weight configuration $\mathbf{W}$ follows the Boltzmann distribution $P(\mathbf{W})\propto e^{-\beta\mathcal{H}(\mathbf{W})}$, where the normalization constant defined by $Z$ is the partition function of the model, and the inverse-temperature $\beta$ determines the energy level one is interested in. We first calculate the partition function of a new system obtained by adding a weight and its associated $C$ constraints (the degree of this weight is $C$), as given by

[TABLE]

where $M$ denotes the total number of active paths (constraints) in the model, $W_{i}$ is the newly added weight, $Z^{old}=\sum_{\mathbf{W}}\exp({\sum_{a=1}^{M}\beta v^{a}\prod_{i\in\partial a}W_{i}})$ is the partition function of the old system, $m_{j\rightarrow b}$ is the cavity magnetization of weight $j$ in the absence of constraint $b$, and $j\in\partial b\backslash i$ denotes the set of weights involved in the constraint $b$ but $i$ is excluded from this set. To derive the second equality in Eq. (5), we use the assumption that when a constraint is removed, its neighboring weights become independent, which is asymptotically correct when the factor graph is locally tree-like. Upon defining the conjugate magnetization $\hat{m}_{b\rightarrow i}\equiv\tanh\beta v^{b}\prod_{j\in\partial b\backslash i}m_{j\rightarrow b}$, one gets the free energy shift:

[TABLE]

Then, we add a constraint (e.g., $a$) to form a new system, and calculate the new partition function as

[TABLE]

The corresponding free energy shift is thus -\beta\Delta F_{a}=\ln\biggl{[}\cosh\beta v^{a}\bigl{(}1+\tanh\beta v^{a}\prod_{i\in\partial a}m_{i\rightarrow a}\bigr{)}\biggr{]}. Finally, the Bethe free energy scaled by $\beta$ can be constructed as follows Mézard and Parisi (2001):

[TABLE]

Following a variational principle of the free energy, one can derive the recursive equation of $\{m_{i\rightarrow a},\hat{m}_{b\rightarrow i}\}$ as follows:

[TABLE]

where $\partial b\backslash i$ denotes the member of interaction $b$ except $i$, and $\partial i\backslash a$ denotes the interaction set $i$ is involved in with $a$ removed. $m_{i\rightarrow a}$ is interpreted as the message passing from the weight $i$ to the interaction $a$ it participates in, while $\hat{m}_{b\rightarrow i}$ is interpreted as the message passing from the interaction $b$ to its member $i$. We call this iteration equation the belief propagation, which serves as the message passing algorithm whose fixed point corresponds to the stationary point of the Bethe free energy, from the variational principle Mézard and Montanari (2009).

With the free energy, one can estimate both the entropy and the energy of the model. The energy is given by:

[TABLE]

The entropy can be obtained by using the standard thermodynamic formula as $S=-F+\beta E$. In the paramagnetic phase, all magnetizations vanish, and analytic expressions for both energy and entropy can be derived. The energy per weight is written as $-\alpha\tanh\beta$, where $\alpha$ denotes the ratio between the number of constructed active paths and the number of weights involved in these paths, and the value of $\alpha$ depends on both $p_{l}$ and $\xi_{l}$, and thus can be estimated after the random construction of the RAP model. Analogously, the entropy per weight is written as $\ln 2+\alpha(\ln\cosh\beta-\beta\tanh\beta)$.

The above equation is derived under the replica symmetry assumption, i.e., the Gibbs state does not split into an exponential number of states where each of them (indexed by $\gamma$) is described by a free energy $F^{\gamma}$. If this assumption is broken, a replica symmetry breaking picture should be introduced. Under this picture, the message passing along each edge (e.g., $i\rightarrow a$) of the factor graph turns out to be a probability density $P_{i\rightarrow a}(m_{i\rightarrow a})$, capturing the fluctuations of a replica symmetric message among different states. The replica symmetry breaking case can be in principle analyzed, yet complicates the analysis, therefore we would not give a full introduction here, and interested readers can find basics in the classic book Mézard and Montanari (2009). We remark here that the replica symmetry solution seems enough to describe the thermodynamic behavior of the investigated model, within the explored range of parameters, which is supported by the stability of the solution and the non-negative entropy for a discrete system.

Appendix B Monte-Carlo method to simulate the RAP model

Given the model defined in Eq. (3), one can explore its energy landscape by using the Monte-Carlo importance sampling method. More precisely, after a random initialization of the weight configuration, the configuration is updated by the following rule: the transition probability from state $\mathbf{W}$ to $\mathbf{W}^{\prime}$ with only $W_{i}$ flipped ($W_{i}^{\prime}=-W_{i}$) is expressed as $\exp(-2\beta W_{i}H_{i})$ where the effective local field is defined as $H_{i}=\sum_{a\in\partial i}v^{a}\prod_{j\in\partial a\backslash i}W_{j}$. In our Monte-Carlo simulations, the system is quenched directly from a very high temperature (random initialization) to the desired temperature. The dynamics evolves for a certain number of sweeps (each sweep consists in $N_{w}$ proposed weight flips), and is finally sampled to evaluate thermodynamic quantities such as energy, spin glass order parameters. We run the dynamics in a total of $10^{4}$ sweeps, and after relaxation, $500$ configurations are collected for numerical computations. The hyper-parameters of Monte-Carlo sampling are chosen to take the trade-off between accuracy and computational cost for a system of $\mathcal{O}(10^{3})\sim\mathcal{O}(10^{4})$ degrees of freedom.

We also perform simulated annealing experiments with $\tau$ Monte-Carlo sweeps per temperature. The temperature is decreased according to the set $\{2T_{f},2T_{f}-\Delta T,\ldots,T_{f}\}$ where $\Delta T=0.1$. We change the value of $\tau$ to study the temperature effects on the energy level that can be reached by the simulated annealing.

Acknowledgements.

We would like to thank Taro Toyoizumi for discussions. We also thank the referee for the stimulating comments. This research was supported by AMED under Grant Number JP18dm020700 and the start-up budget 74130-18831109 of the 100-talent-program of Sun Yat-sen University.