Minimal model of permutation symmetry in unsupervised learning

TL;DR

This paper introduces a minimal two-unit model in a restricted Boltzmann machine to analyze how permutation symmetry influences learning, revealing that correlation between hidden units reduces the data needed for concept formation and that symmetry breaking occurs spontaneously.

Contribution

The study provides a semi-rigorous proof that the critical data size for symmetry breaking is independent of the number of hidden units and introduces an efficient algorithm for inferring hidden unit receptive fields.

Findings

Correlation between hidden units reduces critical data size.

Weakly-correlated receptive fields significantly lower the data threshold.

Permutation symmetry can be spontaneously broken during learning.

Abstract

Permutation of any two hidden units yields invariant properties in typical deep generative neural networks. This permutation symmetry plays an important role in understanding the computation performance of a broad class of neural networks with two or more hidden units. However, a theoretical study of the permutation symmetry is still lacking. Here, we propose a minimal model with only two hidden units in a restricted Boltzmann machine, which aims to address how the permutation symmetry affects the critical learning data size at which the concept-formation (or spontaneous symmetry breaking in physics language) starts, and moreover semi-rigorously prove a conjecture that the critical data size is independent of the number of hidden units once this number is finite. Remarkably, we find that the embedded correlation between two receptive fields of hidden units reduces the critical data size. In particular, the weakly-correlated receptive fields have the benefit of significantly reducing the minimal data size that triggers the transition, given less noisy data. Inspired by the theory, we also propose an efficient fully-distributed algorithm to infer the receptive fields of hidden units. Furthermore, our minimal model reveals that the permutation symmetry can also be spontaneously broken following the spontaneous symmetry breaking. Overall, our results demonstrate that the unsupervised learning is a progressive combination of spontaneous symmetry breaking and permutation symmetry breaking which are both spontaneous processes driven by data streams (observations). All these effects can be analytically probed based on the minimal model, providing theoretical insights towards understanding unsupervised learning in a more general context.

Figures

Equations

$$P(\mathbf{\sigma })=\frac{1}{Z(\mathbf{\xi })}\cosh (\beta X)\cosh (\beta Y),$$

$$\begin{array}{rl}& P(\mathbf{\xi }|\{{\mathbf{\sigma }}^a{\}}_{a=1}^M)=\frac{{\prod }_{a}P({\mathbf{\sigma }}^a|\mathbf{\xi })}{{\sum }_{\mathbf{\xi }}{\prod }_{a}P({\mathbf{\sigma }}^a|\mathbf{\xi })}\\ & =\frac{1}{\Omega }\prod _a\frac{1}{Z({\mathbf{\xi }}^1,{\mathbf{\xi }}^2)}\cosh \left(\frac{\beta }{\sqrt{N}}{\mathbf{\xi }}^1\cdot {\mathbf{\sigma }}^a\right)\cosh \left(\frac{\beta }{\sqrt{N}}{\mathbf{\xi }}^2\cdot {\mathbf{\sigma }}^a\right),\end{array}$$

$$\begin{array}{rl}& Z({\mathbf{\xi }}^1,{\mathbf{\xi }}^2)=\sum _{\mathbf{\sigma }}\cosh \left(\frac{\beta }{\sqrt{N}}{\mathbf{\xi }}^1\cdot \mathbf{\sigma }\right)\cosh \left(\frac{\beta }{\sqrt{N}}{\mathbf{\xi }}^2\cdot \mathbf{\sigma }\right)\\ & =\frac{1}{2}\sum _{\mathbf{\sigma }}\left[\cosh (X+Y)+\cosh (X-Y)\right]\\ & =\frac{1}{2}\left[\prod _{i}2\cosh \left(\frac{\beta }{\sqrt{N}}({\xi }_i^1+{\xi }_i^2)\right)+\prod _{i}2\cosh \left(\frac{\beta }{\sqrt{N}}({\xi }_i^1-{\xi }_i^2)\right)\right]\\ & \simeq 2^N{e}^{{\beta }^2}\cosh ({\beta }^{2}Q),\end{array}$$

$$P(\mathbf{\xi }|\{{\mathbf{\sigma }}^a{\}}_{a=1}^M)=\frac{1}{\Omega }\prod _a\frac{1}{\cosh ({\beta }^{2}Q)}\cosh \left(\frac{\beta }{\sqrt{N}}{\mathbf{\xi }}^1\cdot {\mathbf{\sigma }}^a\right)\cosh \left(\frac{\beta }{\sqrt{N}}{\mathbf{\xi }}^2\cdot {\mathbf{\sigma }}^a\right),$$

$$P_{i\to a}({\xi }_i^1,{\xi }_i^2)$$

$$\begin{array}{rl}{\mu }_{b\to i}({\xi }_i^1,{\xi }_i^2)& =\sum _{\mathbf{\xi }\backslash {\xi }_i^1,{\xi }_i^2}\frac{1}{\cosh \left({\beta }^2{Q}_c+\frac{{\beta }^2}{N}{\xi }_i^1{\xi }_i^2\right)}\cosh \left(\beta X_b+\frac{\beta }{\sqrt{N}}{\xi }_i^1{\sigma }_i^b\right)\cosh \left(\beta Y_b+\frac{\beta }{\sqrt{N}}{\xi }_i^2{\sigma }_i^b\right)\\ & \times \prod _{j\in \partial b\backslash i}{P}_{j\to b}({\xi }_j^1,{\xi }_j^2),\end{array}$$

$$G_{b\to i}^1$$

$$G_{b\to i}^2$$

$${\Gamma }_{b\to i}^1$$

$${\Gamma }_{b\to i}^2$$

$${\Xi }_{b\to i}$$

$$\begin{array}{rl}{\mu }_{b\to i}({\xi }_i^1,{\xi }_i^2)& =\frac{1}{\cosh \left({\beta }^2{Q}_{b\to i}+\frac{{\beta }^2}{N}{\xi }_i^1{\xi }_i^2\right)}\iint DxDy\cosh \left(\beta \sqrt{{\Gamma }_{b\to i}^1}x+\beta G_{b\to i}^1+\frac{\beta }{\sqrt{N}}{\xi }_i^1{\sigma }_i^b\right)\\ & \times \cosh \left(\beta \sqrt{{\Gamma }_{b\to i}^2}(\psi x+\sqrt{1-{\psi }^2}y)+\beta G_{b\to i}^2+\frac{\beta }{\sqrt{N}}{\xi }_i^2{\sigma }_i^b\right),\end{array}$$

$$\begin{split}u_{b\rightarrow i}(\xi_{i}^{1},\xi_{i}^{2})&=\frac{\beta^{2}\Gamma_{b\rightarrow i}^{2}(1-\psi^{2})}{2}-\ln\left(2\cosh\Bigl{(}\beta^{2}Q_{b\rightarrow i}+\frac{\beta^{2}\xi_{i}^{1}\xi_{i}^{2}}{N}\Bigr{)}\right)+\frac{\beta^{2}}{2}\left(\sqrt{\Gamma_{b\rightarrow i}^{1}}+\sqrt{\Gamma_{b\rightarrow i}^{2}}\psi\right)^{2}\\ &+\ln\cosh\left(\beta G_{b\rightarrow i}^{1}+\beta G_{b\rightarrow i}^{2}+\frac{\beta}{\sqrt{N}}\sigma_{i}^{b}(\xi^{1}_{i}+\xi^{2}_{i})\right)\\ &+\ln\left[1+e^{-2\beta^{2}\sqrt{\Gamma_{b\rightarrow i}^{1}\Gamma_{b\rightarrow i}^{2}}\psi}\frac{\cosh\left(\beta G_{b\rightarrow i}^{1}-\beta G_{b\rightarrow i}^{2}+\frac{\beta}{\sqrt{N}}\sigma_{i}^{b}(\xi^{1}_{i}-\xi^{2}_{i})\right)}{\cosh\left(\beta G_{b\rightarrow i}^{1}+\beta G_{b\rightarrow i}^{2}+\frac{\beta}{\sqrt{N}}\sigma_{i}^{b}(\xi^{1}_{i}+\xi^{2}_{i})\right)}\right].\end{split}$$

$$m_{i\to a}^1$$

$$m_{i\to a}^2$$

$$q_{i\to a}$$

$$\Delta f_i$$

$$\begin{array}{rl}\Delta f_a& =\frac{{\beta }^2{\Gamma }_a^2(1-{\tilde{\psi }}^2)}{2}-\ln \left(2\cosh ({\beta }^2{Q}_a)\right)+\frac{{\beta }^2}{2}{\left(\sqrt{{\Gamma }_a^1}+\sqrt{{\Gamma }_a^2}\tilde{\psi }\right)}^2\\ & +\ln \cosh \left(\beta G_a^1+\beta G_a^2\right)+\ln \left[1+e^{-2{\beta }^2{\Xi }_a}\frac{\cosh \left(\beta G_a^1-\beta G_a^2\right)}{\cosh \left(\beta G_a^1+\beta G_a^2\right)}\right],\end{array}$$

$$-\beta f=\underset{n\to 0,N\to \infty }{\lim }\frac{\ln ⟨{\Omega }^n⟩}{nN},$$

$$\begin{array}{rl}⟨{\Omega }^n⟩& =\sum _{\{{\mathbf{\sigma }}^a,{\mathbf{\xi }}^{\mathrm{true}}\}}\prod _i\left[P({\xi }_i^{1,\mathrm{true}},{\xi }_i^{2,\mathrm{true}})\right]\prod _a\frac{\cosh \left(\frac{\beta }{\sqrt{N}}{\mathbf{\xi }}^{1,\mathrm{true}}\cdot {\mathbf{\sigma }}^a\right)\cosh \left(\frac{\beta }{\sqrt{N}}{\mathbf{\xi }}^{2,\mathrm{true}}\cdot {\mathbf{\sigma }}^a\right)}{2^N{e}^{{\beta }^2}\cosh ({\beta }^{2}q)}\\ & \times \sum _{\{{\mathbf{\xi }}^{1,\gamma },{\mathbf{\xi }}^{2,\gamma }\}}\prod _{a,\gamma }\frac{\cosh \left(\frac{\beta }{\sqrt{N}}{\mathbf{\xi }}^{1,\gamma }\cdot {\mathbf{\sigma }}^a\right)\cosh \left(\frac{\beta }{\sqrt{N}}{\mathbf{\xi }}^{2,\gamma }\cdot {\mathbf{\sigma }}^a\right)}{\cosh ({\beta }^2{R}^{\gamma })},\end{array}$$

$$\begin{split}-\beta f_{{\rm RS}}&=-R\hat{R}-T_{1}\hat{T}_{1}-\tau_{1}\hat{\tau}_{1}-T_{2}\hat{T}_{2}-\tau_{2}\hat{\tau}_{2}+\frac{\hat{q}_{1}(q_{1}-1)}{2}+\frac{\hat{q}_{2}(q_{2}-1)}{2}+\frac{r\hat{r}}{2}\\ &+\alpha\beta^{2}\left(1-\frac{q_{1}+q_{2}}{2}\right)-\alpha\ln\Bigl{(}2\cosh(\beta^{2}R)\Bigr{)}+\int D\mathbf{z}\left[\ln Z_{{\rm eff}}\right]_{\xi^{1,{\rm true}},\xi^{2,{\rm true}}}\\ &+\frac{\alpha e^{-\beta^{2}}}{\cosh(\beta^{2}q)}\int D\mathbf{t}\cosh(\beta t_{0})\cosh(\beta qt_{0}+\beta\sqrt{1-q^{2}}x_{0})\ln Z_{{\rm E}},\end{split}$$

$$Z_{\mathrm{eff}}$$

$$b_1$$

$$b_2$$

$$b_3$$

$$Z_{\mathrm{E}}$$

$${\Lambda }_{+}$$

$${\Lambda }_{-}$$

$$T_1$$

$$T_2$$

$$q_1$$

$$q_2$$

$${\tau }_1$$

$${\tau }_2$$

$$R$$

$$r$$

$${\hat{T}}_1$$

$${\hat{T}}_2$$

$${\hat{q}}_1$$

$${\hat{q}}_2$$

$${\hat{\tau }}_1$$

$${\hat{\tau }}_2$$

$$\hat{R}$$

$$\hat{r}$$

Full text

1 Introduction

Unsupervised learning is defined as the process of searching for latent features in raw (unlabeled) data, and thus serves as a fundamental property of the cerebral cortex of the brain [1, 2]. To understand unsupervised learning from a neural network perspective, restricted Boltzmann machine (RBM) is proposed. RBM is a two-layered neural network with one layer called the visible layer, and the other called the hidden layer. No lateral connections exist in each layer. The connections between visible and hidden layers are called synaptic weights, representing encoded latent features in the observed data. The process of learning the synaptic weights from the unlabeled data (also called training) mimics the unsupervised learning. RBM is thus receiving substantial research interests both from machine learning and statistical physics communities [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13].

Training of RBM relies on the maximum likelihood principle via a gradient ascent procedure. The mean activity of each neuron or correlations between visible and hidden neurons can be estimated by either truncated Gibbs sampling [3, 14] or advanced mean-field methods [6, 7]. However, the gradient ascent method is difficult to analyze and thus not amenable for a theoretical model. Therefore, based on the probabilistic graphical model framework, one-bit RBM where only one hidden neuron is considered was proposed to address a fundamental issue of unsupervised learning [15], i.e., how many data samples are needed for a successful learning. This work revealed a continuous spontaneous symmetry breaking (SSB) transition separating a random-guess phase from a concept-formation phase at a critical value of the amount of provided samples (data size) [16], which is similar to the retarded learning phase transition observed in a generalized Hopfield model of pattern learning [17]. This conclusion is later generalized to RBM with generic priors [18, 19], and synapses of ternary values [20]. However, it is still challenging to handle the case of multiple hidden neurons from the perspective of understanding the learning process as a phase transition. In the presence of multiple hidden neurons, permutation symmetry appears, i.e., the model of the observed data is invariant with respect to exchange of arbitrary two hidden neurons. In addition, the permutation symmetry is a common feature in many modern neural network architectures [21]. Therefore, understanding how the permutation symmetry affects the concept-formation process is important, which may provide us core mechanisms of unsupervised learning.

Here, we propose a minimal model of the permutation symmetry in unsupervised learning, based on mean-field approximations. We show that it is possible to theoretically understand the permutation symmetry using physics approximations. To be more precise, we consider a RBM with two hidden neurons, and embed a latent feature that generates a certain number of data samples through Gibbs samplings of the original model [15]. Then, the data samples are learned by a theory-inspired algorithm, and finally the learned synaptic weights (a latent feature vector) are compared with the embedded ones, to test whether SSB applies to the minimal model, and in addition investigate key factors affecting the critical data size for learning and moreover how the permutation symmetry affects the learning process. We first apply the cavity approximation in statistical mechanics of disordered systems [22] to derive the learning algorithm from a Bayesian inference perspective, whose computation performances in single instances of the model are then predicted by a replica theory. This theory introduces many copies of the original model, and the interaction between any two copies is characterized by a set of self-consistent mean-field equations, from which the critical data size for learning is determined, and moreover whether the permutation symmetry can be spontaneously broken is clarified.

2 Model definition and mean-field methods

2.1 Minimal model of permutation symmetry

In this study, we use the RBM defined above with two hidden neurons (Fig. 1) to learn embedded features in input data samples, which are raw unlabeled data. Each data sample is specified by an Ising-like spin configuration $\bm{\sigma}=\{\sigma_{i}=\pm 1\}_{i=1}^{N}$ where $N$ is the input dimensionality. A collection of $M$ samples is denoted as $\{\bm{\sigma}^{a}\}_{a=1}^{M}$. Synaptic values connecting visible and hidden neurons are characterized by $\bm{\xi}$, where each component takes a binary value ($\pm 1$) as well. Because of two hidden neurons, $\bm{\xi}=(\bm{\xi}^{1},\bm{\xi}^{2})$ where the superscript indicates the hidden neuron’s index. $\bm{\xi}^{1}$ and $\bm{\xi}^{2}$ are also called receptive fields of the first and second hidden neurons, respectively. Statistical properties of this RBM are thus described by the Boltzmann distribution [6]

[TABLE]

where $X=\frac{1}{\sqrt{N}}\bm{\xi}^{1}\cdot\bm{\sigma}$, $Y=\frac{1}{\sqrt{N}}\bm{\xi}^{2}\cdot\bm{\sigma}$, and $Z(\bm{\xi})$ is the partition function depending on the feature $\bm{\xi}$. Note that the two hidden neurons’ activities ($\pm 1$) have been marginalized out. The scaling factor $\frac{1}{\sqrt{N}}$ ensures that the argument of the hyperbolic cosine function is of the order of unity. $\bm{\sigma}$ can be arbitrary one of the $M$ samples. When the embedded feature is randomly generated, the inverse-temperature $\beta$ tunes the noise level of generated data samples from the feature. Clearly, the data distribution is invariant with respect to (w.r.t) the exchange of the hidden neurons, which is called the permutation symmetry in this paper. The required number of hidden neurons to yield this symmetry is at least two, therefore, this setup defines a minimal model to study the permutation symmetry in unsupervised learning.

In this model, the embedded feature follows the distribution $P(\bm{\xi})=P(\bm{\xi}^{1})P(\bm{\xi}^{2}|\bm{\xi}^{1})$ in which $P(\bm{\xi}^{1})=\prod^{N}_{i=1}\left[\frac{1}{2}\delta(\xi_{i}^{1}-1)+\frac{1}{2}\delta(\xi^{1}_{i}+1)\right]$ together with $P(\bm{\xi}^{2}|\bm{\xi}^{1})=\prod_{i=1}^{N}\left[p_{\rm d}\delta(\xi^{2}_{i}=-\xi_{i}^{1})+(1-p_{\rm d})\delta(\xi_{i}^{2}=\xi_{i}^{1})\right]$, where $p_{\rm d}$ controls the fraction of components taking different values in the two feature maps associated with the two hidden neurons.

Given the $M$ data samples, one gets the posterior probability of the embedded feature according to the Bayes’ rule:

[TABLE]

where $\Omega$ is the partition function of the minimal model. For simplicity, a uniform prior for $\bm{\xi}$ is assumed, i.e., we have no prior knowledge about $\bm{\xi}$, although there may exist correlations between two feature maps. In addition, we use the same temperature as that used to generate data. Because we do not use the true prior $\prod_{i}P_{i}(\xi^{1}_{i},\xi^{2}_{i}|p_{{\rm d}})$, the current setting does not require the value of $p_{{\rm d}}$, thereby is not the Bayes-optimal setting which corresponds to Nishimori condition in physics [16]. Therefore, using the uniform prior is more computationally challenging. We leave a detailed analysis of the Bayes-optimal setting in a future work.

One obstacle to compute the posterior probability is the nested partition function $Z(\bm{\xi}^{1},\bm{\xi}^{2})$. Fortunately, this partition function can be simplified in the large-$N$ limit. More precisely,

[TABLE]

where we have used $\ln\cosh(x)\simeq\frac{x^{2}}{2}$ for small $x$ to arrive at the final equality, and $Q\equiv\frac{1}{N}\sum_{i}\xi_{i}^{1}\xi_{i}^{2}$, which is exactly the overlap between the two feature maps. To sum up, we move all the irrelevant constants into the partition function $\Omega$, the posterior probability can then be rewritten as

[TABLE]

which forms the Boltzmann distribution of our minimal model. In this paper, we consider the case of $M=\alpha N$ where $\alpha$ specifies the data (constraint) density.

2.2 Cavity approximation to handle the posterior probability

In what follows, we compute the maximizer of the posterior marginals (MPM) estimator $(\hat{\xi}^{1}_{i},\hat{\xi}_{i}^{2})=\arg\max_{\xi^{1}_{i},\xi_{i}^{2}}P_{i}(\xi_{i}^{1},\xi_{i}^{2})$ [16], where the feature map of each hidden neuron is combined and the prediction is thus the augmented version of the inferred feature vector in the one-bit RBM [15]. Hence, the task is to compute marginal probabilities, i.e., $P_{i}(\xi_{i}^{1},\xi_{i}^{2})$, which is intractable due to the interaction among data constraints (the product over $a$ in Eq. (4)). However, by mapping the original model (Eq. (4)) onto a graphical model (Fig. 1), where data constraints and paired-synapses are treated respectively as factor (data) nodes and variable nodes, one can estimate the marginal probability by running a message passing iteration among factor and variable nodes, as we shall explain below. The key assumption is that the paired-synapses on the graphical model are weakly correlated, which is called the Bethe approximation [23] in physics.

We first define a cavity probability $P_{i\rightarrow a}(\xi_{i}^{1},\xi_{i}^{2})$ with the data node $a$ removed. Under the weak correlation assumption, $P_{i\rightarrow a}(\xi_{i}^{1},\xi_{i}^{2})$ obeys a self-consistent equation:

[TABLE]

where $Z_{i\rightarrow a}$ is a normalization constant, $\partial i\backslash a$ denotes neighbors of the feature node $i$ except the data node $a$, $\partial b\backslash i$ denotes neighbors of the data node $b$ except the feature node $i$, and the auxiliary quantity $\mu_{b\rightarrow i}(\xi_{i}^{1},\xi_{i}^{2})$ denotes the contribution from data node $b$ to feature node $i$ given the value of $(\xi_{i}^{1},\xi_{i}^{2})$ [6, 15]. Products in Eq. (5) result from the weak correlation assumption. In addition, $X_{b}\equiv\frac{1}{\sqrt{N}}\sum_{j\neq i}\xi^{1}_{j}\sigma_{j}^{b}$, $Y_{b}\equiv\frac{1}{\sqrt{N}}\sum_{j\neq i}\xi^{2}_{j}\sigma_{j}^{b}$, and the cavity version of $Q$ is defined as $Q_{c}\equiv\frac{1}{N}\sum_{j\neq i}\xi_{j}^{1}\xi_{j}^{2}$, which can be further replaced by its typical value obtained by the average over the cavity probability (to be shown below). Although this is a crude approximation, it works quite well in practice.

Still, the above self-consistent equation is intractable due to the summation to estimate $\mu_{b\rightarrow i}$. Nevertheless, a careful inspection reveals that $X_{b}$ and $Y_{b}$ are approximately correlated Gaussian random variables due to the central limit theorem. As a result, the intractable summation can be replaced by an integral which is easy to calculate in this model. We just need to compute the following mean, variance and covariance between these random variables.

[TABLE]

where $G$ and $\Gamma$ denotes the mean and variance of the Gaussian random variable respectively, and the last quantity denotes the covariance between $X_{b}$ and $Y_{b}$. The cavity magnetization is defined as $m_{j\rightarrow b}^{1,2}=\sum_{\xi_{j}^{1},\xi_{j}^{2}}\xi_{j}^{1,2}P_{j\rightarrow b}(\xi_{j}^{1},\xi_{j}^{2})$, and the cavity correlation is defined as $q_{j\rightarrow b}=\sum_{\xi_{j}^{1},\xi_{j}^{2}}\xi_{j}^{1}\xi_{j}^{2}P_{j\rightarrow b}(\xi_{j}^{1},\xi_{j}^{2})$. Finally, using the above parameters of the correlated Gaussian distribution, we rewrite $\mu_{b\rightarrow i}(\xi_{i}^{1},\xi_{i}^{2})$ as

[TABLE]

where $Dx\equiv\frac{e^{-x^{2}/2}dx}{\sqrt{2\pi}}$, $\psi=\frac{\Xi_{b\rightarrow i}}{\sqrt{\Gamma_{b\rightarrow i}^{1}\Gamma_{b\rightarrow i}^{2}}}$, and $Q_{b\rightarrow i}=\frac{1}{N}\sum_{j\neq i}q_{j\rightarrow b}$ stemming from $Q_{c}$ in Eq. (5b) replaced by its cavity mean. The above integral representation of $\mu_{b\rightarrow i}(\xi^{1}_{i},\xi^{2}_{i})$ can be analytically estimated; for convenience, we define $u_{b\rightarrow i}(\xi_{i}^{1},\xi_{i}^{2})\equiv\ln\mu_{b\rightarrow i}(\xi_{i}^{1},\xi^{2}_{i})$. It is easy to show that

[TABLE]

To close the iteration equation, we need to compute the cavity magnetization and correlation as follows:

[TABLE]

$m_{i\rightarrow a}^{1,2}$ can be interpreted as the message passing from feature node $i$ to data node $a$ ($q_{i\rightarrow a}$ is also similarly interpreted), while $u_{b\rightarrow i}$ can be interpreted as the message passing from data node $b$ to feature node $i$.

If the weak correlation assumption is self-consistent, starting from randomly initialized messages, the learning equations will converge to a fixed point corresponding to a thermodynamically dominant minimum of the Bethe free energy function [24], which is given by $-\beta f_{{\rm Bethe}}=\frac{1}{N}\sum_{i}\Delta f_{i}-\frac{N-1}{N}\sum_{a}\Delta f_{a}$. The free energy contributions of variable node and data node are given respectively by:

[TABLE]

where $\tilde{\psi}=\frac{\Xi_{a}}{\sqrt{\Gamma_{a}^{1}\Gamma_{a}^{2}}}$. The forms of $\Gamma_{a}^{1,2}$, $G_{a}^{1,2}$, $Q_{a}$ and $\Xi_{a}$ are similar to their cavity counterparts (e.g., in Eq. (8)), but with the only difference that the node $i$’s contribution is not excluded. Once the iteration converges, the MPM estimator predicts that $\hat{\xi}^{1}_{i}={\rm sgn}(m_{i}^{1})$ and $\hat{\xi}^{2}_{i}={\rm sgn}(m_{i}^{2})$, where the full (non-cavity) magnetization $m_{i}^{1,2}$ is computed taking into account all contributions of adjacent data nodes to the node $i$ (see Eq. (9), and the symbol $\backslash a$ is thus removed).

2.3 Replica theory of the minimal model

To have an analytic argument about the critical threshold for spontaneous symmetry breaking, we calculate the free energy in the thermodynamic limit using the replica method. Instead of calculating a disorder average of $\ln\Omega$, the replica method computes the disorder average of an integer power of $\Omega$, i.e., $\left<\Omega^{n}\right>$. In physics, this corresponds to preparing $n$ replicas of the original system; then the rescaled free energy density (multiplied by $-\beta$) can be obtained as [16]

[TABLE]

where the limits of $N\rightarrow\infty$ and $n\rightarrow 0$ have been exchanged, such that the thermodynamic limit can be taken first for applying the Laplace’s method or saddle-point analysis [25], and the disorder average is taken over all possible samplings (data) and the random realizations of the true feature vector. The explicit form of $\left<\Omega^{n}\right>$ reads

[TABLE]

where $\gamma$ indicates the replica index, $\bm{\xi}^{{\rm true}}\equiv(\bm{\xi}^{1,{\rm true}},\bm{\xi}^{2,{\rm true}})$, $q\equiv\frac{1}{N}\bm{\xi}^{1,{\rm true}}\cdot\bm{\xi}^{2,{\rm true}}$, and $R^{\gamma}\equiv\frac{1}{N}\bm{\xi}^{1,\gamma}\cdot\bm{\xi}^{2,\gamma}$. Note that $q$ is pre-determined and used to generate the random true feature maps, as also defined in section 2.1. We leave the technical details to A, and give the final result here. The free energy function reads,

[TABLE]

where $[\cdot]_{\xi^{1,{\rm true}},\xi^{2,{\rm true}}}$ means an average w.r.t $P(\xi^{1,{\rm true}},\xi^{2,{\rm true}})$, $D\mathbf{z}\equiv Dz_{1}Dz_{2}Dz_{3}$ (a standard Gaussian measure vector, as defined below Eq. (7)), and similarly $D\mathbf{t}\equiv Dt_{0}Dx_{0}DuDu^{\prime}$. ${\rm RS}$ means the replica symmetry assumption we used to get the final result. This assumption implies that the order parameter (various kinds of overlaps, explicitly defined below) does not rely on its specific replica index. We assume that this assumption is able to describe the system as we shall show it leads to consistent predictions verified in algorithmic results of single instances. The auxiliary quantities $Z_{{\rm eff}}$ and $Z_{{\rm E}}$ are defined as follows,

[TABLE]

where $A\equiv T_{1}\tau_{1}+\frac{(\tau_{2}-T_{1}q)(T_{2}-\tau_{1}q)}{1-q^{2}}$, $B\equiv\sqrt{q_{1}-T_{1}^{2}-\frac{(\tau_{2}-T_{1}q)^{2}}{1-q^{2}}}$, and $K\equiv\sqrt{q_{2}-\tau_{1}^{2}-\frac{(T_{2}-\tau_{1}q)^{2}}{1-q^{2}}-\frac{(r-A)^{2}}{B^{2}}}$.

The associated (non-conjugated) saddle-point equations are expressed as

[TABLE]

Note that the average w.r.t the true features can be written explicitly by definition as $P(\xi^{1,{\rm true}},\xi^{2,{\rm true}})=\frac{p_{{\rm d}}}{2}$ for both true components taking different values, and otherwise $P(\xi^{1,{\rm true}},\xi^{2,{\rm true}})=\frac{1-p_{{\rm d}}}{2}$. $p_{{\rm d}}$ is related to $q$ by $p_{{\rm d}}=\frac{1-q}{2}$. The outer average also includes the disorder average over $\mathbf{z}$. The inner average $\left<\bullet\right>$ indicates the thermal average under the partition function $Z_{{\rm eff}}$ (corresponding to a two-spin interaction Hamiltonian). This average is analytically tractable, e.g., $\left<\xi^{1}\right>=\frac{1}{Z_{{\rm eff}}}\frac{\partial Z_{{\rm eff}}}{\partial b_{1}}=\frac{\tanh b_{1}+\tanh b_{2}\tanh b_{3}}{1+\tanh b_{1}\tanh b_{2}\tanh b_{3}}$. $\left<\xi^{2}\right>$ and $\left<\xi^{1}\xi^{2}\right>$ can also be similarly computed.

We further comment that $T_{1}$ characterizes the typical overlap between inferred value of the first feature and its true counterpart, and likewise $T_{2}$ characterizes the typical overlap between the second feature and its ground truth; $q_{1}$ and $q_{2}$ characterize the sizes of the first and second feature spaces respectively; $R$ characterizes the correlation of the two features within the same replica, while $r$ is the correlation for different replicas; $\tau_{1}$ characterizes the typical correlation between the first true feature and the inferred value of the second feature in an arbitrary replica, and likewise $\tau_{2}$ characterizes the typical correlation between the second true feature and the inferred value of the first feature in an arbitrary replica. $\tau_{1}$ and $\tau_{2}$ are thus responsible for the permutation symmetry effect. Taken all together, $(T_{1},T_{2},q_{1},q_{2},R,r,\tau_{1},\tau_{2})$ forms the order parameter set of our model. Their exact mathematical definitions are given in the A.

Finally, the conjugated order parameters can also be derived from a saddle point analysis of the free energy function, and they obey the following equations:

[TABLE]

where the average $\left<\bullet\right>\equiv\frac{e^{-\beta^{2}}}{\cosh\beta^{2}q}\int D\mathbf{t}\cosh(\beta t_{0})\cosh(\beta qt_{0}+\beta\sqrt{1-q^{2}}x_{0})\bullet$, $\langle\langle\bullet\rangle\rangle\equiv\frac{e^{-\beta^{2}}}{\cosh\beta^{2}q}\int D\mathbf{t}\sinh(\beta t_{0})\cosh(\beta qt_{0}+\beta\sqrt{1-q^{2}}x_{0})\bullet$, and $\langle\langle\langle\bullet\rangle\rangle\rangle\equiv\frac{e^{-\beta^{2}}}{\cosh\beta^{2}q}\int D\mathbf{t}\cosh(\beta t_{0})\sinh(\beta qt_{0}+\beta\sqrt{1-q^{2}}x_{0})\bullet$. The auxiliary quantities are defined as follows,

[TABLE]

To sum up, Eqs. (15) and (16) construct a closed iterative equation (detailed derivations are given in the B), whose fixed point gives an approximate evaluation of the free energy. When all order parameters vanish (a trivial disordered state), the free energy has an analytic value expressed as $\alpha\beta^{2}+\ln 4$ in agreement with $-\beta f_{\rm Bethe}$ in the same trivial state. In addition, these saddle point equations in the case of $q=0$ can be simplified to the result of Ref. [16] for unsupervised feature learning in a one-bit RBM (see the C). More precisely, when the true feature maps are orthogonal, we have $\Omega=\Omega_{{\rm one-bit-RBM}}^{2}$, thus the free energy is two times as large as that of one-bit RBM.

Meanwhile, the converged order parameters from Eqs. (15) and (16), especially $T_{1}$ and $T_{2}$ can be compared with the algorithmic results, and can also be used to analytically derive the critical threshold $\alpha_{c}$ for unsupervised learning in this permutation-symmetry model. When the data size is not sufficient, we expect that the order parameters vanish, and in the small order-parameter limit, $T_{1}\simeq\hat{T}_{1}+q\hat{\tau}_{2}$, $\tau_{2}\simeq\hat{\tau}_{2}+\hat{T}_{1}q$, $\hat{T}_{1}\simeq\alpha\beta^{4}[T_{1}+\tau_{2}\tanh(\beta^{2}q)]$, and $\hat{\tau}_{2}\simeq\alpha\beta^{4}[T_{1}\tanh(\beta^{2}q)+\tau_{2}]$. Based on these four equations, it is easy to show that the critical learning threshold is given by

[TABLE]

As expected, the critical value for our current model does not depend on the sign of $q$. A detailed derivation of $\alpha_{c}$ is given in the D. In the correlation-free limit $q=0$, the known result of the one-bit RBM, $\alpha_{c}=\beta^{-4}$ [15, 16] is recovered. Thus, we first theoretically prove the conjecture made by empirical observations in Ref. [18] that once the number of hidden neurons is finite (not proportional to $N$), the critical learning threshold does not change with this finite number! However, Ref. [18] overlooks the effect of the potential correlation across true hidden features, which indeed affects the learning threshold. Remarkably, this effect, more natural than the ideal correlation-free case, is clearly captured by our theory (Eq. (18)). We show this effect in Fig. 2. We can see that for a fixed noise level, increasing $q$ has the effect of decreasing the critical data size. That is, if the data set is created by strongly correlated receptive fields (feature maps), then the learning is relatively easy, or fewer samples are sufficient to trigger the phase transition of concept-formation. Conversely, if $q$ is zero, the data is generated from independent feature maps, then a successful learning requires a much more larger dataset. We also observe that for a large $\beta$ (less noisy data), a small correlation level $q$ can already significantly reduce the minimal data size that triggers learning, as verified by the fact that around the small $q$ region the larger the $\beta$ is, the sharper the surface becomes.

Next, we analyze two interesting limits implied by the critical threshold equation (Eq. (18)). In the limit $|q|\rightarrow 1$, $\alpha_{c}\rightarrow\frac{1}{4}\beta^{-4}$ provided that $\beta$ is relatively large such that $\tanh\beta^{2}\simeq 1$. The second case is another limit $|q|\rightarrow 0$, i.e., $q$ takes a small value but not zero, implying that a weak correlation among feature maps is maintained. Depending on the order of magnitude of $q$, we have the following result given a relatively large $\beta$:

[TABLE]

Note that $\infty$ means any large value of $\beta$ such that $\tanh\beta\simeq 1$ rather than a definite value of infinity. Eq. (19) reveals that once the two feature maps are weakly correlated, the minimal learning data size for a transition can be further (or even significantly) reduced compared to the correlation-free case, particularly in the case that $q$ is not very small but still larger than the order of magnitude set by $\beta^{-2}$. We show this result in Fig. 3.

We thereby have a significant hypothesis for the triggering of concept-formation that a bit large (compared with $\beta^{-2}$) yet still small value of the correlation level is highly favored for unsupervised learning from a dataset of smaller size (compared with the correlation-free case). Regularization techniques such as locally enforcing feature orthogonality [26] and dropping some connections during training [27] have been introduced to deep learning. Weakly-correlated receptive fields are also favored from the perspective of neural computation, since the redundancy among synaptic weights is reduced and thus different feature detectors inside the network can encode efficiently stimuli features rather than capturing noise in the data. A similar decorrelation in hidden activities was recently theoretically analyzed in feedforward neural networks [28]. We hope our theoretical prediction can be verified in specific machine learning tasks, and even in neuroscience experiments where the relationship amongst the minimal data size for learning, the correlation level of synapses (or receptive fields) and the noise level in stimuli can be jointly established. Therefore, from the Bayesian learning perspective, the correlated-feature-map case yields a much lower threshold of phase transition towards the concept formation, in comparison with the correlation-free case  [15, 16, 18].

Overall, the prediction quantitatively captures the learning behavior in both uncorrelated and correlated settings. In next section, we shall further verify this conclusion by extensive numerical simulations on single instances of the minimal model.

3 Results and discussion

In this section, we study how the permutation symmetry between two hidden neurons affects the learning process, i.e., the spontaneous symmetry breaking transition of the concept-formation during unsupervised learning. We focus on whether the replica theory predicts the learning threshold related to the continuous transition. The learning threshold can be estimated from the message passing algorithmic results on single instances of the minimal model. We first randomly generate true feature maps with a pre-determined correlation level specified by $q$. The true feature maps are then used to generate $M$ Monte-Carlo samples through a Gibbs sampling procedure [15] according to Eq. (1). Finally, these samples are used as a quenched disorder for the Bayesian inference of the true feature maps [16]. An overlap with the ground truth is computed and compared with the replica prediction. In fact, for comparison, we define the overlap [16], e.g., $T_{1}^{\rm{MP}}=\frac{1}{N}\sum_{i}\xi_{i}^{\rm{true}}m_{i}^{1}$ where $\rm{MP}$ means message passing, and $m_{i}^{1}$ takes into account the thermal average (uncertainty about the ground truth). Other overlaps can be similarly defined.

First, we compare the free energy function estimated under the Bethe approximation with that predicted by the replica theory (in the thermodynamic limit). As shown in Fig. 4 (a), the algorithmic results on finite-sized networks coincide very well with the theoretical predictions. This implies that, the approximation we used to derive the message passing equation for learning a RBM with two hidden units is reasonable, especially when the number of visible neurons is large.

Then, we study the evolution of the overlap with the ground truth as a function of the data density. As shown in Fig. 4 (b), we clearly see a continuous phase transition separating a disordered (symmetric) phase from an ordered (symmetry-broken) phase. This is consistent with previous works revealing the spontaneous symmetry breaking in unsupervised feature learning [15, 16]. When the amount of provided data samples is not too large, the original model maintains its symmetry (i.e., equal probabilities for positive and negative assignments of synapses respectively). This is in an isotropic phase which does not capture any concept from the data samples. However, the increasing amount of data samples will break this symmetry through a continuous phase transition towards a non-trivial concept formation. During this process, in a practical message passing procedure, messages flowing in the factor graph are biased towards the true feature maps underlying the noisy data. Remarkably, the theory predicts the exact location of the phase transition point at which the messages running on a single instance start to polarize towards the true feature map. As predicted and observed, the learning threshold indeed decreases as the absolute value of $q$ grows. For $|q|=0.3$, $\alpha_{c}\simeq 0.596$; and for $|q|=0.8$, $\alpha_{c}$ is significantly reduced to be about $0.334$.

Due to the permutation symmetry, we find that after the SSB transition, there appear in sequence three non-trivial solutions to the saddle point equations (Eqs. (15) and (16)), as $\alpha$ increases for a given value of $(\beta,q)$. The first solution has the form of $q_{1}=q_{2}=r$ and $T_{1}=T_{2}=\tau_{1}=\tau_{2}$. This solution is caused by the permutation symmetry and is dominant at the earlier stage of the post-SSB-transition. $q_{1}=q_{2}=r$ implies the permutation symmetry between inferred vectors of $\bm{\xi}^{1}$ and $\bm{\xi}^{2}$; and $T_{1}=T_{2}=\tau_{1}=\tau_{2}$ implies the permutation symmetry between planted true feature vectors. In other words, after the SSB transition, the unsupervised learning increases the overlap with the ground truth through identifying the common components of the two true feature maps.

The SSB phase is stable until a point where the unsupervised learning starts to predict the different components (as shown in Fig. 5 (a)), thereby breaking the symmetry between the inferred values of $\bm{\xi}^{1}$ and $\bm{\xi}^{2}$. This point is thus referred to as $\alpha_{c}^{{\rm PSB_{s}}}$, namely, the permutation symmetry breaking (PSB) for student/inferred features. However, the inferred features are equally likely to be either $(\bm{\xi}^{1,{\rm true}},\bm{\xi}^{2,{\rm true}})$ or its permutation, thus the permutation symmetry between planted/teacher features (i.e., $\bm{\xi}^{1,{\rm true}}$ and $\bm{\xi}^{2,{\rm true}}$) is still preserved, therefore, after the ${\rm PSB_{s}}$ transition, we have the second solution: $q_{1}=q_{2}\neq r$, and $T_{1}=T_{2}=\tau_{1}=\tau_{2}$. This second solution is stable until there is a turnover of the trend of the permutation-type overlap or the order parameter $r$ (as shown in the insets of Fig. 4 (b) and Fig. 5 (a) respectively). At this turnover, $r$ starts to decrease, thereby reducing the permutation-type overlap towards the true correlation level of the two feature maps (Fig. 4 (b)). After this transition indicated by the turnover, even the permutation symmetry between the teacher’s features is broken, as the unsupervised learning is able to distinguish the two feature maps underlying the raw data samples. Therefore, we observe the final solution that has the form of either $q_{1}=q_{2}=T_{1}=T_{2}$ or $q_{1}=q_{2}=\tau_{1}=\tau_{2}$ (two sub-forms), which is dominant at the later stage of the unsupervised learning. Accordingly, this second PSB transition is called ${\rm PSB_{t}}$ transition with broken symmetry for teacher’s features.

We remark that this third type non-trivial solution can be deduced from the fact that as $\alpha$ grows, the inferred feature map gets close to the true feature map and thus both the inferred and true feature maps follow the same posterior probability of the learning process. Note that for this solution, these two sub-forms share the same free energy. One inferred feature map has the freedom of choosing to match the first or second true feature map. This choice does not change the overall free energy. Therefore, at the later stage, one solution of a larger overlap corresponds to the case that the inferred feature maps are matched with their true counterparts; the other smaller overlap corresponds to the case that the inferred feature maps are matched with their interchanged (permuted) counterparts. In particular, the latter one should converge to the embedded correlation $q$ as the learning converges to the true feature maps. We show one example in the case of $q=0.3$ (see the inset of Fig. 4 (b)), which is consistent with the replica prediction (the inset of Fig. 5 (a)). Note that these two kinds of overlaps are very hard to distinguish at the earlier stage, since the inferred feature maps are only partially correct compared with their true counterparts. Indeed, in the algorithmic simulations, we observe two solutions of different overlaps (one is larger than the other), which are very easy to distinguish at the later stage of the learning. In addition, considering the sign of each inferred (or true) feature vector, we have other multiple types of solutions, since a sign-reversed assignment of the feature map does not change the posterior probability (Eq. (4)). For simplicity, we do not consider the sign-reversed case in the analysis.

Permutation symmetry of different types can be spontaneously broken. Both transitions are continuous. We summarize these qualitatively different transitions in Fig. 5 (b), showing that as $q$ grows, the gap between $\alpha_{c}^{{\rm SSB}}$ and $\alpha_{c}^{{\rm PSB}}$ increases, thereby demonstrating that the small-$q$ region is beneficial not only in the sense of a significant decrease of data size triggering the concept-formation but also in the sense of a small data size triggering the permutation symmetry breaking. Therefore, not only can increasing the number of data samples (observations) drive a spontaneous symmetry breaking, but increasing the observations also drive a permutation symmetry breaking which finally leads to a perfect reconstruction of the embedded feature maps.

In addition, when one tries to run the message passing from a perturbed true feature map, it is observed that its basin is quite large, i.e., it is stable in the presence of a sufficiently large dataset. But if the perturbation is large enough (e.g., crossing the sign boundary), the passing message is able to converge to other types of fixed points (e.g., the sign-reversed or permutation-symmetric ones).

4 Summary

In conclusion, we propose a minimal model of the permutation symmetry in a simple unsupervised learning system— a simple RBM with only two hidden units. Using statistical mechanics tools developed in theory of disordered systems, we reveal very rich properties of this model. Effects of permutation symmetry among hidden units have been studied in supervised learning in multilayer neural networks [29, 30]. Here we are the first to consider its effects on unsupervised learning, focusing on spontaneous symmetry breaking and the critical data size that triggers the transition.

First, we derive an efficient algorithm to infer the embedded feature vectors from a given dataset (Monte-Carlo samples in this paper), according to the Bayesian principle. This algorithm is fully distributed, in that the computation is implemented in terms of local passing messages among feature nodes and data (constraint) nodes. The algorithmic results can be used as a test of the replica theory in the thermodynamic limit.

Second, the behavior of the algorithm, especially the critical data size at which a continuous phase transition from a disordered (symmetric) phase to an ordered (symmetry-broken) phase occurs, can be predicted by a replica theory of the minimal model. Due to the existence of two hidden units, we have to manipulate eight order parameters that characterize all possible typical correlations in the replica space itself as well as between replica space and true-feature space. In addition to the spontaneous symmetry breaking, our theory predicts another later transition where the permutation symmetry is spontaneously broken and thus the identity of the feature map can be captured during unsupervised learning. By studying the minimal model, we are able to interpret the unsupervised learning, a fundamental process governing artificial and biological intelligence, as a progressive combination of SSB and PSB (including two types—${\rm PSB_{s}}$ and ${\rm PSB_{t}}$), both of which are driven by increasing the data size (or observations).

Based on the replica computation, we have three contributions. (i) We analytically prove that the critical data size for the phase transition indeed does not depend on the number of hidden units (a finite number case), and this conclusion was empirically observed in a previous work [18]. The critical value $\alpha_{c}=\beta^{-4}$ for learning in a one-bit RBM [15, 16] can directly apply to the RBM with two hidden units, once no correlations are embedded into the two feature vectors of hidden units. A detailed proof is also given in the C. (ii) In addition, we reveal the correlation level of embedded true feature vectors reduces the critical data size, characterized by a simple formula (see Eq. (18)). As an example, in the very large $\beta$ limit (small variability in the data space), the critical data size for $q=-1$ is reduced to only one-fourth of that in the correlation-free case. In another limit $q\rightarrow 0$, depending on the order of magnitude of $q$ compared with $\beta^{-2}$, the necessary data size that triggers the transition can be reduced to one half of $\beta^{-4}$ or even less, although the receptive fields are weakly-correlated. This prediction qualitatively coincides with the observation that humans or non-human animals do not need a large amount of data samples to learn a concept from structured examples (e.g, natural images are highly structured as various levels of correlations are embedded) [31]. We expect that this quantitative prediction of the critical learning data threshold can be shown to hold in a generic case with arbitrary levels of receptive-field-correlation, and with an arbitrary finite number of hidden units, in particular for generative models of neural networks. (iii) Our theory predicts that an additional spontaneous permutation symmetry breaking follows the spontaneous symmetry breaking, leading to another benefit of the small-$q$ regime where the data size triggering the permutation symmetry breaking is smaller compared with the large-$q$ regime (Fig. 5).

Our study also encourages several interesting future directions. First, by analogy with a Bayesian iteration derived in our previous work [16] to predict the true noise level (the hyper-parameter $\beta$) in a dataset, one can also derive the iteration equation for predicting both correlation level $q$ and noise level $\beta$ for an arbitrary dataset. Second, one can also verify the predicted value of $\alpha_{c}$ in a practical neural network architecture with inferred $q$ and $\beta$. This can be carried out in a more complex RBM, by comparing the predicted value of the critical data size with the observed one for a successful unsupervised learning task. Furthermore, the prediction of our theory about the benefit of weakly-correlated feature maps compared with the correlation-free situation can also be tested in practical artificial neural networks and even in biological neural networks. Finally, it would be very interesting to verify in a more general setting whether the progressive combination of SSB and PSB is the underlying mechanism of unsupervised learning. In this sense, our minimal model paves the way towards understanding the fundamentally important unsupervised learning process, by addressing the role of the permutation symmetry and moreover the critical learning data sizes (related to both SSB and PSB) in a simple unsupervised learning system.

Acknowledgments

We would like to thank referees for their inspiring comments. This research was supported by the start-up budget 74130-18831109 of the 100-talent- program of Sun Yat-sen University (H.H.), the NSFC (Grant No. 11805284) (H.H.), and grants from the Research Grants Council of Hong Kong (16322616 and 16396817) (T.H. and M.W.)

Appendix A Computation of $\left<\Omega^{n}\right>$

Here, we show details to compute $\left<\Omega^{n}\right>$, which is defined as:

[TABLE]

where $\gamma$ indicates the replica index, $\bm{\xi}^{true}=\{\bm{\xi}^{1,true},\bm{\xi}^{2,true}\}$, the overlap $q=\frac{1}{N}\bm{\xi}^{1,true}\bm{\xi}^{2,true}$, and $R^{\gamma}=\frac{1}{N}\bm{\xi}^{1,\gamma}\bm{\xi}^{2,\gamma}$.

To further calculate $\langle\Omega^{n}\rangle$, we need to define the order parameters as follows:

[TABLE]

Note that these order parameters can be used to evaluate the disorder average in Eq. (20). We then proceed as follows,

[TABLE]

where we have defined $X^{0}_{a}=\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\xi_{i}^{1,true}\sigma_{i}^{a}$, $Y^{0}_{a}=\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\xi_{i}^{2,true}\sigma_{i}^{a}$, and $X_{a}^{\gamma}=\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\xi_{i}^{1,\gamma}\sigma_{i}^{a}$, $Y^{\gamma}_{a}=\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\xi_{i}^{2,\gamma}\sigma_{i}^{a}$. To get the second equality, we have used the integral representation of the delta function $\delta(x)=\int\frac{d\hat{x}}{2\pi}e^{i\hat{x}x}$.

To compute the free energy value, we assume a simple ansatz, i.e., all order parameters do not depend on their specific replica indexes, which is called the replica-symmetry assumption. To be more precise, we assume

[TABLE]

for any $\gamma$. We also assume that

[TABLE]

for any $\gamma$ and $\gamma^{{}^{\prime}}$. Then we can express $\langle\Omega^{n}\rangle$ as :

[TABLE]

In the thermodynamics limit, $\langle\Omega^{n}\rangle$ can be approximated as $e^{N\mathcal{A}(\mathcal{O}_{*},\hat{\mathcal{O}}_{*},\alpha,\beta,n)}$(namely, the saddle-point method), where $\mathcal{O}_{*}$ and $\mathcal{\hat{O}}_{*}$ represent all non-conjugated order parameters and conjugated order parameters evaluated at the maximal value of the action, respectively. The expression for the action $\mathcal{A}(\mathcal{O},\hat{\mathcal{O}},\alpha,\beta,n)$ can be written by

[TABLE]

where $G_{S}$ is the entropy term, and $G_{E}$ is the energy term.

To derive the entropy term $G_{S}$, we use the following identities[32]:

[TABLE]

The above non-linear terms can be reduced to linear terms in the exponential functions of Eq. (25) by the Hubbard-Stratonovich transformation $\int Dte^{bt}=e^{\frac{1}{2}b^{2}}$. Then, we obtain $G_{S}$ as :

[TABLE]

Finally, we can express the entropy term $G_{S}$ in a compact form as

[TABLE]

where we have defined $D\mathbf{z}=Dz_{1}Dz_{2}Dz_{3}$, and the auxiliary variables $b_{1},b_{2},$ and $b_{3}$ as

[TABLE]

We remark that in the expression of $G_{S}$, the inner summation over $\xi^{1},\xi^{2}$ can be thought as a two-spin interaction partition function, which is defined as $Z_{{\rm eff}}$ in the main text. $[\bullet]_{\xi^{1,true},\xi^{2,true}}$ means an average w.r.t $P(\xi^{1,true},\xi^{2,true})$ which is also defined in the main text.

Next, we turn to compute the energy term $G_{E}$. The expression of $G_{E}$ is given by

[TABLE]

where $\langle\bullet\rangle$ defines the disorder average. $X^{0}$, $Y^{0}$, $X^{\gamma}$, $Y^{\gamma}$ are correlated Gaussian random variables, which are the same as before but data index $a$ has been dropped off. They have zero mean and unit variance. Their covariances are determined by the aforementioned order parameters as follows:

[TABLE]

The random variables $X^{0},Y^{0},X^{\gamma},Y^{\gamma}$ can thus be parameterized by six standard Gaussian variables of zero mean and unit variance ($t_{0},x_{0},u,u^{{}^{\prime}},y_{\gamma},\omega_{\gamma}$) as follows,

[TABLE]

where $A=T_{1}\tau_{1}+\frac{(\tau_{2}-T_{1}q)(T_{2}-\tau_{1}q)}{1-q^{2}}$, $B=\sqrt{q_{1}-(T_{1})^{2}-\frac{(\tau_{2}-T_{1}q)^{2}}{1-q^{2}}}$, and $K=\sqrt{q_{2}-(\tau_{1})^{2}-\frac{(T_{2}-\tau_{1}q)^{2}}{1-q^{2}}-(\frac{r-A}{B})^{2}}$. One can easily verify that the above parameterization satisfies their covariance structures. Therefore, the $G_{E}$ term can be calculated by a standard Gaussian integration given by

[TABLE]

where $C\equiv\sqrt{1-q_{2}-\frac{(R-r)^{2}}{1-q_{1}}}$. To proceed, we first define the auxiliary quantities as :

[TABLE]

Then we compute the integral inside the power $n$ defined by $I$ whose result is given by

[TABLE]

For simplicity, we also define the following auxiliary quantities $Z_{E},G_{c}^{-},G_{s}^{+},G_{s}^{-}$:

[TABLE]

Following the replica trick, we can get:

[TABLE]

where the integral in the denominator can be exactly computed with the result given by

[TABLE]

Finally, by collecting all the above relevant terms, we have the following estimation of $\left<\Omega^{n}\right>$ given by

[TABLE]

where in shorthand $D\mathbf{t}=Dt_{0}Dx_{0}DuDu^{\prime}$. By computing $\lim_{n\to 0}\frac{\ln{\langle\Omega^{n}\rangle}}{n}$ and using Eq. (38), we get the expression $F_{\beta}=-\beta f_{RS}$ as

[TABLE]

Note that we have used $\lim_{n\to 0}\frac{\ln\left[{\int D\mathbf{z}Z^{n}_{{\rm eff}}}\right]_{\xi^{1,true},\xi^{2,true}}}{n}=\int D\mathbf{z}\left[\ln Z_{{\rm eff}}\right]_{\xi^{1,true},\xi^{2,true}}$ to arrive at the final expression.

Appendix B Derivation of saddle-point equations

By the saddle-point analysis, these non-conjugated order parameters $\mathcal{O}$ should obey the following stationary conditions:

[TABLE]

Similarly, for conjugated order parameters $\mathcal{\hat{O}}$, the following stationary conditions should be satisfied as:

[TABLE]

We first evaluate the self-consistent equations those non-conjugated order-parameters obey. For $R$, we have the following equation as

[TABLE]

Thus the saddle-point equation of $R$ is given by

[TABLE]

where the thermal average $\langle\bullet\rangle$ is computed under the partition function $Z_{{\rm eff}}$ (a two-spin interaction partition function), and the outer average indicates the disorder average over Gaussian random variables $\mathbf{z}$ and the distribution $P(\xi^{1,true},\xi^{2,true})$.

Similarly, for the order parameter $T_{1}$, we have the following equation as

[TABLE]

Noting that $\frac{\partial Z_{eff}}{\partial\hat{T_{1}}}=\sum_{\xi^{1},\xi^{2}}\xi^{1,true}\xi^{1}e^{b_{1}\xi^{1}+b_{2}\xi^{2}+b_{3}\xi^{1}\xi^{2}}$, we get the the final expression of $T_{1}$ as

[TABLE]

The expressions of $T_{2},\tau_{1}$ and $\tau_{2}$ can be derived in the same way as follows:

[TABLE]

[TABLE]

[TABLE]

Next, we turn to the saddle-point equation of $q_{1}$, i.e.,

[TABLE]

Noticing that $\frac{\partial Z_{eff}}{\partial\hat{q}_{1}}=\frac{1}{2}(\hat{q_{1}}-\frac{\hat{r}}{2})^{-\frac{1}{2}}\sum_{\xi^{1},\xi^{2}}\xi^{1}z_{1}e^{b_{1}\xi^{1}+b_{2}\xi^{2}+b_{3}\xi^{1}\xi^{2}}$, we get the expression of $q_{1}$ as

[TABLE]

To proceed, we use the following identity

[TABLE]

where $f(z)$ is any differentiable function of $z$. Thus we have the following equality as

[TABLE]

Finally, the expression of $q_{1}$ is given by

[TABLE]

Similarly, $q_{2}$ should obey the following equation given by

[TABLE]

Following the same spirit, we get the following stationary condition for $r$ as

[TABLE]

Note that

[TABLE]

By applying Eq. (53), we can obtain the following three identities as

[TABLE]

Using the above three identities together with Eq. (54), we get the expression of the saddle-point equation for $r$ as follows

[TABLE]

Given the result that $Z_{eff}=2e^{b_{3}}\cosh{(b_{1}+b_{2})}+2e^{-b_{3}}\cosh{(b_{1}-b_{2})}$, the thermal average like $\langle\xi^{1}\rangle,\langle\xi^{2}\rangle$, and $\langle\xi^{1}\xi^{2}\rangle$ can be easily calculated as follows:

[TABLE]

and

[TABLE]

and finally

[TABLE]

In case of $\hat{r}<0$, we can re-parameterize $b_{1}$ and $b_{2}$ as

[TABLE]

We remark that this re-parameterization does not change the final results of multidimensional Gaussian integrations in the saddle-point equation.

To sum up, the saddle-point equations for non-conjugated order parameters are given by

[TABLE]

Next, we derive the saddle-point equations for those conjugated order parameters. For $\hat{R}$, we obtain the saddle point equation as

[TABLE]

where $\frac{\partial}{\partial R}\ln{Z_{E}}=\beta^{2}\frac{e^{\beta^{2}(R-r)}\cosh{(\beta\Lambda_{+})}-e^{-\beta^{2}(R-r)}\cosh{(\beta\Lambda_{-}})}{e^{\beta^{2}(R-r)}\cosh{(\beta\Lambda_{+})}+e^{-\beta^{2}(R-r)}\cosh{(\beta\Lambda_{-})}}=\beta^{2}G_{c}^{-}.$ Therefore, the saddle-point equation of $\hat{R}$ is given by

[TABLE]

For convenience, we define the measure $\langle\bullet\rangle$ as $\frac{e^{-\beta^{2}}}{\cosh{(\beta^{2}}q)}\int D\mathbf{t}\cosh{\beta t_{0}}\cosh{\beta(qt_{0}+\sqrt{1-q^{2}}x_{0})}\bullet$. As a result,

[TABLE]

For $\hat{T_{1}}$, we have the following condition

[TABLE]

To proceed, we first get the derivation of $\Lambda_{+}$ and $\Lambda_{-}$ w.r.t $T_{1}$ as follows

[TABLE]

Then, the derivation of $\ln{Z_{E}}$ w.r.t $T_{1}$ can be simplified into the form as

[TABLE]

To further simplify the result, we need to evaluate the following equations. The first one is derived by applying Eq. (53) as

[TABLE]

The second one is derived as

[TABLE]

The third one is derived as

[TABLE]

The fourth one is derived as

[TABLE]

The last one is given by

[TABLE]

Through a bit lengthy algebraic manipulations, we get from Eq. (69)

[TABLE]

We thus define another measure $\langle\langle\bullet\rangle\rangle=\frac{e^{-\beta^{2}}}{\cosh{(\beta^{2}q})}\int D\mathbf{t}\sinh{\beta t_{0}}\cosh{\beta(qt_{0}+\sqrt{1-q^{2}}x_{0})}\bullet$, and it then follows that

[TABLE]

Similarly, we can obtain the saddle-point equation of $\hat{\tau_{1}}$ as

[TABLE]

Next we turn to the saddle-point equations for $\hat{T}_{2}$ and $\hat{\tau_{2}}$. We first get the derivation of $\Lambda_{+}$ and $\Lambda_{-}$ w.r.t $T_{2}$ as

[TABLE]

Based on the above equations, we get the derivation of $\ln{Z_{E}}$ w.r.t $T_{2}$ given by

[TABLE]

Then we have

[TABLE]

For a further simplification, we need to derive the following identity as

[TABLE]

Using Eq. (83) together with Eq. (75) and Eq. (76), we finally arrive at the saddle-point equation of $\hat{T}_{2}$:

[TABLE]

We thus define the third measure $\langle\langle\langle\bullet\rangle\rangle\rangle=\frac{e^{-\beta^{2}}}{\cosh{(\beta^{2}q})}\int D\mathbf{t}\cosh{\beta t_{0}}\sinh{\beta(qt_{0}+\sqrt{1-q^{2}}x_{0})}\bullet$. We then write the saddle-point equation in a compact form as

[TABLE]

Similarly, we obtain the saddle-point equation for $\hat{\tau_{2}}$ as

[TABLE]

Then we turn to the saddle-point equations of $\hat{q}_{1}$ and $\hat{q}_{2}$. From $\frac{\partial F_{\beta}}{\partial q_{1}}=0$, we get

[TABLE]

The derivation of $\ln{Z_{E}}$ w.r.t $q_{1}$ is given by

[TABLE]

Using Eq. (74), Eq. (75) and Eq. (76), we get the saddle-point equation of $\hat{q}_{1}$ as

[TABLE]

Similarly, we can derive the saddle-point equation for $\hat{q}_{2}$ as

[TABLE]

Lastly, we derive the saddle-point equation for $\hat{r}$ as

[TABLE]

Noting that $\frac{\partial\ln{Z_{E}}}{\partial r}=-\beta^{2}G_{c}^{-}+\beta\left(\frac{1}{B}G_{s}^{-}u+\frac{\partial K}{\partial r}G_{s}^{-}u^{{}^{\prime}}\right)$, we get the saddle-point equation of $\hat{r}$ as

[TABLE]

Note that Eq. (75) and Eq. (76) are used to derive the final result.

To sum up, the saddle-point equations of our minimal model are listed as follows

[TABLE]

Appendix C The free energy function in the limit of $q=0$

In the case of $q=0$, the saddle point equation of the minimal model has the solution: $q_{1}=q_{2}=T_{1}=T_{2}$ and other order parameters vanish. Thus, we can simplify $\Lambda_{+}$ and $\Lambda_{-}$ as follows,

[TABLE]

We then define $\chi_{1}=T_{1}t_{0}+\sqrt{q_{1}-(T_{1})^{2}}u$, and $\chi_{2}=T_{2}x_{0}+\sqrt{q_{2}-(T_{2})^{2}}u^{\prime}$, so the saddle point eqtation of $\hat{T}_{1}$ is given by

[TABLE]

where we used the identity $\int Dx_{0}\cosh(\beta x_{0})=e^{\beta^{2}/2}$. Similarly, one can prove that $\hat{T}_{1}=\hat{T}_{2}$. As for $\hat{q}_{1}$, we will have

[TABLE]

Similarly, one can prove that $\hat{q}_{1}=\hat{q}_{2}$.

It is easy to prove that $\hat{\tau}_{1}=0,\hat{\tau_{2}}=0$ and $\hat{R}=0,\hat{r}=0$, then we can express $b_{1},b_{2}$ and $b_{3}$ as :

[TABLE]

Therefore, $\int D\mathbf{z}[\ln Z_{eff}]_{\xi^{1,true},\xi^{2,true}}$ can be simplified as $2\int Dz\ln 2\cosh(\hat{T}_{1}+\sqrt{\hat{q}_{1}}z)$. $T_{1}$ becomes

[TABLE]

One can easily prove that $T_{1}=T_{2}$. Similarly for the order parameter $q_{2}$, we can also get:

[TABLE]

It is easy to show that $q_{1}=q_{2}$, and moreover $R=r=\tau_{1}=\tau_{2}=0$. To sum up, we recover the saddle point equations of one-bit RBM reported in Ref. [16].

Next, we show the $q=0$ version of the free energy function. It is easy to show that $Z_{{\rm E}}=\cosh{\beta(\chi_{1}+\chi_{2})}+\cosh{\beta(\chi_{1}-\chi_{2})}=2\cosh{\beta\chi_{1}}\cosh{\beta\chi_{2}}$. Therefore, we have the following integral

[TABLE]

Collecting all the relevant terms, we can show that the free energy of our minimal model with $q=0$ is two times as large as that of one-bit RBM, which can also be intuitively understood by the argument that the partition function factorizes as $\Omega=\Omega_{{\rm one-bit-RBM}}^{2}$. Therefore we can conclude that the critical data size for spontaneous symmetry breaking does not change even if an additional hidden node is added. This conclusion seems to carry over to the case of more hidden nodes following the principle of the partition function’s factorization.

Appendix D Derivation of the critical data size (Eq. (18) in the main text)

We assume that near to the transition point, all order parameters are very small such that we can expand them to leading order. According to Eq. (65), when the critical point is approached from below, $\langle\xi^{1}\rangle\simeq\tanh{b_{1}}\simeq b_{1}$. Analogously, $\langle\xi^{2}\rangle\simeq b_{2}$, and $\langle\xi^{1}\xi^{2}\rangle\simeq b_{3}$. We thus have the following equalities in this limit:

[TABLE]

Similarly, in the limit of vanishing order parameters, we have the following approximation

[TABLE]

Inserting this approximation into the saddle-point equations of $\hat{T}_{1}$ and $\hat{\tau}_{2}$, we obtain the approximate results of $\hat{T_{1}}$ and $\hat{\tau_{2}}$ as

[TABLE]

We recast the equations for all these four order parameters in a matrix form as

[TABLE]

From the Eq. (111) and Eq. (118), $T_{1}$ and $\tau_{2}$ can be calculated out as

[TABLE]

where the matrix $\mathcal{M}$ is named the stability matrix, whose largest eigenvalue determines the critical value of the learning data size $\alpha_{c}$. In detail, the stability matrix has two eigenvalues:

[TABLE]

The $\alpha_{c}$ can be read off from $\lambda_{+}=1$, i.e.,

[TABLE]

An alternative way to understand that the smaller eigenvalue could not be used to determine $\alpha_{c}$, is that it leads to a non-physical solution $\alpha_{c}=\frac{\beta^{-4}}{1+q\tanh{\beta^{2}q}-|q+\tanh{(\beta^{2}q)}|}$. Because in a special case of large $\beta$ limit and positive $q$, $\tanh{(\beta^{2}q)}\simeq 1-2e^{-2\beta^{2}q}$, then we have $\alpha_{c}\simeq\frac{e^{2\beta^{2}q}}{2(1-q)\beta^{4}}$, which implies that this value tends to $\infty$ which is in contradiction with the expectation that learning should be easier given noise-free data.