Energy--Information Trade-off Induces Continuous and Discontinuous Phase Transitions in Lateral Predictive Coding

TL;DR

This paper analyzes how energy and information trade-offs in lateral predictive coding lead to phase transitions in neural network structures, revealing insights into the emergence of complex internal models.

Contribution

It provides an analytical and numerical study of phase transitions in lateral predictive coding, highlighting the effects of energy-information trade-offs on network structure and efficiency.

Findings

Identification of continuous and discontinuous phase transitions in synaptic weights

Discovery of a reciprocity-breaking transition building cyclic dominance

Optimal networks saturate energy efficiency bounds

Abstract

Lateral predictive coding is a recurrent neural network which creates energy-efficient internal representations by exploiting statistical regularity in sensory inputs. Here we investigate the trade-off between information robustness and energy in a linear model of lateral predictive coding analytically and by numerical minimization of a free energy. We observe several phase transitions in the synaptic weight matrix, especially a continuous transition which breaks reciprocity and permutation symmetry and builds cyclic dominance and a discontinuous transition with the associated sudden emergence of tight balance between excitatory and inhibitory interactions. The optimal network follows an ideal-gas law in an extended temperature range and saturates the efficiency upper-bound of energy utilization. These results bring theoretical insights on the emergence and evolution of complex internal models in predictive processing systems.

Equations

$$\frac{\mathrm{d}{x}_i}{\mathrm{d}t}=s_i-x_i-\sum _{j\ne i}{w}_{ij}{x}_j,$$

$$E\equiv\sum_{\vec{\bm{s}}}{\rm P}(\vec{\bm{s}}){\vec{\bm{x}}}_{\vec{\bm{s}}}^{2}=\textrm{Tr}\Bigl{[}\bigl{(}\bm{I}+\bm{W}\bigr{)}^{-1}\bm{C}\bigl{(}\bm{I}+\bm{W}^{\top}\bigr{)}^{-1}\Bigr{]}\;,$$

$$S=-\log\bigl{[}{\rm Det}(\bm{I}+\bm{W})\bigr{]}$$

$$F=E-TS,$$

$$T_{N}^{\textrm{rb}}=2\bigl{(}\sqrt{1+(N-1)c}+(N-1)\sqrt{1-c}\bigr{)}^{2}/N^{2}\;.$$

$$E=\frac{N}{2}T(N\ge 3)$$

$$\mathbf{W}=\left[\begin{array}{ccccc}0& w_1& w_2& \cdots & w_{N-1}\\ w_{N-1}& 0& w_1& \cdots & w_{N-2}\\ w_{N-2}& w_{N-1}& 0& \cdots & w_{N-3}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ w_1& w_2& w_3& \cdots & 0\end{array}\right].$$

$$\mathcal{O}_{CD}=1-\min\limits_{(p_{1},\ p_{2},\ p_{N})}\biggl{|}\frac{\sum_{i=1}^{N-1}w_{p_{i+1}\ p_{i}}}{\sum_{i=1}^{N-1}w_{p_{i}\ p_{i+1}}}\biggr{|}\;,$$

$$\mathcal{O}_{EI}=\frac{1}{N}\sum\limits_{i=1}^{N}\biggl{[}1-\frac{\bigl{|}\sum_{j\neq i}w_{ij}\bigr{|}}{\sum_{j\neq i}|w_{ij}|}\biggr{]}\;.$$

$$x_i=g_i(\mathbf{s})+{\xi }_i,$$

$$g_{i}(\bm{s})=\sum_{j=1}^{N}\Bigl{(}\frac{\bm{I}}{\bm{I}+\bm{W}}\Bigr{)}_{ij}s_{j}\;,$$

$$\mathcal{I}[\mathbf{x};\mathbf{s}]=H[\mathbf{x}]+H[\mathbf{s}]-H[\mathbf{x}|\mathbf{s}],$$

$$P_{\text{output}}(\mathbf{x})$$

$$\approx$$

$$H[\bm{x}]\equiv-\int{\rm d}\bm{x}P_{\textrm{output}}(\bm{x})\ln P_{\textrm{output}}(\bm{x})=-\ln\Bigl{|}\textrm{Det}\bigl{(}\bm{I}+\bm{W}\bigr{)}\Bigr{|}+H[\bm{s}]\;.$$

$$x_1=\frac{s_1-w_{12}{s}_2}{1-w_{12}{w}_{21}},x_2=\frac{s_2-w_{21}{s}_1}{1-w_{12}{w}_{21}}.$$

$$E=\bigl{\langle}(x_{1})^{2}\bigr{\rangle}+\bigl{\langle}(x_{2})^{2}\bigr{\rangle}=\frac{c_{11}+c_{22}-2c_{12}\,(w_{12}+w_{21})+c_{22}\,w_{12}^{2}+c_{11}\,w_{21}^{2}}{(1-w_{12}\,w_{21})^{2}}\;.$$

$$c_{11}=⟨s_1^2⟩,c_{22}=⟨s_2^2⟩,c_{12}=⟨s_1{s}_2⟩.$$

$$S=-\ln\bigl{(}1-w_{12}\,w_{12}\bigr{)}\;,$$

$$\mathbf{C}=(1-c)\mathbf{I}+c\left[\begin{array}{c}1\\ 1\\ ⋮\\ 1\end{array}\right]\left[\begin{array}{cccc}1& 1& \dots & 1\end{array}\right],$$

$$S=-\ln\bigl{(}1+(N-1)w\bigr{)}-(N-1)\ln(1-w)\;.$$

$$\bm{M}=\frac{\bm{I}}{\bm{I}+\bm{W}}\bm{C}\Bigl{(}\frac{\bm{I}}{\bm{I}+\bm{W}}\Bigr{)}^{\top}\;.$$

$$\bm{M}=\frac{1-c}{(1-w)^{2}}\,\bm{I}+\frac{1}{N}\Bigl{(}\frac{1+(N-1)c}{(1+(N-1)w)^{2}}-\frac{1-c}{(1-w)^{2}}\Bigr{)}\,\begin{bmatrix}1\\ 1\\ \vdots\\ 1\end{bmatrix}\,\begin{bmatrix}1\;&1\;&\ldots\;&1\end{bmatrix}\;,$$

$$E\equiv\textrm{Tr}\bigl{[}\bm{M}\bigr{]}=\frac{1+(N-1)c}{\bigl{(}1+(N-1)w\bigr{)}^{2}}+(N-1)\frac{1-c}{(1-w)^{2}}\;.$$

$$T=\frac{2}{Nw}\Bigl{(}\frac{(1-c)\bigl{(}1+(N-1)w\bigr{)}}{(1-w)^{2}}-\frac{(1-w)\bigl{(}1+(N-1)c\bigr{)}}{[1+(N-1)w]^{2}}\Bigr{)}\;.$$

$$\frac{\mathbf{I}}{\mathbf{I}+\mathbf{W}+\mathbf{\epsilon }}=\frac{\mathbf{I}}{\mathbf{I}+\mathbf{W}}-\frac{\mathbf{I}}{\mathbf{I}+\mathbf{W}}\mathbf{\epsilon }\frac{\mathbf{I}}{\mathbf{I}+\mathbf{W}}+\frac{\mathbf{I}}{\mathbf{I}+\mathbf{W}}\mathbf{\epsilon }\frac{\mathbf{I}}{\mathbf{I}+\mathbf{W}}\mathbf{\epsilon }\frac{\mathbf{I}}{\mathbf{I}+\mathbf{W}}+\dots .$$

$$\Delta E$$

$$\frac{\partial S}{\partial w_{ij}}$$

$$\frac{{\partial }^{2}S}{\partial w_{ij}\partial w_{kl}}$$

$$\Delta S=-\textrm{Tr}\Bigl{[}\bm{\varepsilon}\frac{\bm{I}}{\bm{I}+\bm{W}}\Bigr{]}+\frac{1}{2}\textrm{Tr}\Bigl{[}\bm{\varepsilon}\frac{\bm{I}}{\bm{I}+\bm{W}}\bm{\varepsilon}\frac{\bm{I}}{\bm{I}+\bm{W}}\Bigr{]}\;.$$

$$\Delta F^{(1)}=-\textrm{Tr}\Bigl{[}\bm{M}\bm{\varepsilon}^{\top}\frac{\bm{I}}{\bm{I}+\bm{W}^{\top}}+\frac{\bm{I}}{\bm{I}+\bm{W}}\bm{\varepsilon}\bm{M}\Bigr{]}+T\,\textrm{Tr}\Bigl{[}\bm{\varepsilon}\frac{\bm{I}}{\bm{I}+\bm{W}}\Bigr{]}\;.$$

Full text

††preprint: ITP307/2022-03

.1 Introduction

Predictive coding is an influential theory in computational neuroscience with a long history tracing back to Kant and Helmhotz [1, 2, 3, 4], and its recent conceptual advances include the Bayesian brain theory and the free energy principle of information processing and decision making [5, 6, 7, 8]. The predictive coding theory regards the brain as a multi-layered hierarchical neuron network which builds an internal model for the external world and employs it to perceive sensory signals and make predictions. So far most studies on predictive coding have focused on the feedforward and feedback interactions between different layers, and have largely ignored lateral (horizontal, recurrent) interactions within the same layer of neurons. But lateral interactions are abundant in biological neural networks and they play significant roles in perception and inference [9, 10, 11, 12]. Children brain consumes $50\%$ of the body’s total metabolic energy and this fraction decreases to $20\%$ in adulthood [13]. As a big energy consumer it is vital for the brain to be capable of achieving energy-efficient information processing [14, 15, 16, 17, 18]. Lateral predictive coding (LPC) might have been an elegant evolutionary solution to cope with the brain’s tight energy budget.

One basic function of lateral interactions is to reduce the energy cost of internal representations. By adapting the synaptic weights according to the input correlation between spatially adjacent neurons, the external drive to one neuron is partially cancelled by the messages from its neighboring neurons and its activity becomes weaker and sparser [19, 20, 21, 22]. On the other hand, reducing the magnitude and redundancy of internal representation compromises its information content, making it less robust to external and internal noises [23, 24, 25]. How to balance the conflicting demands of energy reduction and information robustness is an important issue of LPC. Yet LPC has been a rarely discussed topic in the physics community, in stark contrast to the huge and enduring theoretical enthusiasm on the other basic neural network models such as the perceptron, the Hopfield model of associative memory and the restricted Boltzmann machine [26].

Here we study phase transitions induced by energy–information trade-off in the simplest model of linear LPC [27]. We introduce a temperature parameter $T$ to incorporate the fitness effect of information robustness and define a free energy, which binds together energy and function, as the minimization objective [28, 29, 30]. We find that even if the input signals are statistically symmetric among all the $N$ units, the optimal synaptic weight matrix will spontaneously break permutation (reciprocal) symmetry at certain critical temperature and form cyclic-dominant interaction patterns. An ideal-gas law $E=(N/2)T$ is then followed by the energy $E$ of the multi-unit ($N\geq 3$) LPC networks, achieving upper-bound efficiency in energy utilization. At low temperatures the weight matrix experiences several further phase transitions, one of them being a discontinuous transition to tight balance of excitatory and inhibitory interactions. Internal symmetries within sub-groups of units may also break, leading to multi-level functional differentiation.

Our theoretical results reveal the emergence of structural complexity in LPC systems. These collective properties of non-reciprocity, correlation reduction, and excitation–inhibition (EI) balance are qualitatively similar to the empirical observations on real-life recurrent neural networks. As an initial attempt we have left many important issues untouched, such as nonlinear synaptic interactions and non-quadratic energetic cost. More efforts from the physics community are needed to fully appreciate LPC as a basic neural computing circuit.

.2 Model and free energy

Consider a fully connected network of $N$ units (each may be a single neuron or a small region of the brain). Denote an internal state vector as $\vec{\bm{x}}=(x_{1},\ldots,x_{N})^{\top}$ and a sensory input as $\vec{\bm{s}}=(s_{1},\ldots,s_{N})^{\top}$. Upon receiving an input $\vec{\bm{s}}$, unit $i$ responds by changing its state $x_{i}$ according to a linear LPC dynamics [19, 21, 27],

[TABLE]

where positive (negative) synaptic weight $w_{ij}$ means that unit $j$ inhibits (excites) unit $i$. The steady-state encoding from sensory input to internal representation is $\vec{\bm{x}}_{\vec{\bm{s}}}=(\bm{I}+\bm{W})^{-1}\vec{\bm{s}}$ with $\bm{I}$ being the identity matrix and $\bm{W}$ the synaptic weight matrix. The convergence of the dynamics (1) requires the real part of every eigenvalue of $(\bm{I}+\bm{W})$ be positive (we set a threshold value $10^{-5}$ and call it the eigenvalue bottomline). We interpret $\sum_{j\neq i}w_{ij}x_{j}$ as the predicted sensory input to unit $i$, and $x_{i}$ as the prediction error between the actual $s_{i}$ and the predicted value [9]. Assuming the energy cost of maintaining and transmitting a prediction error be quadratic in $x_{i}$, the mean energy $E$ of an internal state is

[TABLE]

where ${\rm P}(\vec{\bm{s}})$ is the probability distribution of sensory inputs and $\bm{C}$ is the input correlation matrix with elements $c_{ij}\equiv\sum_{\vec{\bm{s}}}{\rm P}(\vec{\bm{s}})s_{i}s_{j}$. The total energy cost of representing $\mathcal{M}\gg 1$ sensory inputs is then $\mathcal{M}E$.

An infinitesimal volume of the sensory space is transformed by $\vec{\bm{x}}_{\vec{\bm{s}}}$ to an infinitesimal volume in the internal representation space, with Jacobian $1/{\rm Det}(\bm{I}+\bm{W})$ which is the inverse determinant. The entropy of the internal states $\bm{x}$ is then

[TABLE]

plus a constant, and the mutual information between $\bm{s}$ and $\bm{x}$ is equal to this entropy up to a constant [31]. We can take $\mathcal{M}S$ as the total entropy of $\mathcal{M}$ internal states. It is desirable for the entropy $S$ to be large so that the encoding will be sensitive to variations in $\vec{\bm{s}}$ and be robust to noises in representing and transmitting $\vec{\bm{x}}$ (the information maximization principle [28]).

It turns out that minimizing energy (cost) $E$ and maximizing entropy (information) $S$ are mutually conflicting goals. We introduce a temperature parameter $T$ to quantify this energy–information trade-off. In the thermodynamic limit of infinite number of sensory inputs ($\mathcal{M}\rightarrow\infty$), the optimization objective is the total free energy $\mathcal{M}F$ with the coefficient $F$ being

[TABLE]

which combines both energetic and entropic effects [28, 29, 30]. The temperature embodies all the external and internal fitness stresses. At high $T$ values information sensitivity and robustness is the main fitness drive and the adaptation of synaptic weights will favor entropy $S$; at the other limit of low $T$, reducing energy $E$ becomes the dominant fitness concern and the weight matrix will evolve towards energy minimization. The thermodynamic equation accompanying free-energy minimization is simply $T={{\rm d}E}/{{\rm d}S}$.

We search for the global minimum of $F$ at each fixed value of $T$ by a slow simulated annealing process and determine the corresponding optimal weight matrix [31]. We find that $E(T)$ and $S(T)$ are singular at several critical points of $T$, and the optimal weight matrix $\bm{W}$ changes its qualitative properties at these points. Discontinuities of the energy susceptibility (${\rm d}E/{\rm d}T$) imply continuous phase transitions, and discontinuities in the energy itself manifest discontinuous phase transitions. To demonstrate most explicitly the endogenous nature of such phase transitions, we will focus on symmetric and homogeneous sensory inputs: the self-correlation $c_{ii}=1$ for all the inputs $s_{i}$ and pair correlation $c_{ij}=c$ being the same ($c<1$) for the inputs of any two different units. The default internal model $\bm{W}$ is then reciprocal ($w_{ij}=w_{ji}$) and permutation-symmetric (identical non-diagonal elements) [31].

The phase diagram for the simplest two-unit system (Fig. 1) is a good starting point to understand optimal LPC. There is a continuous reciprocity-breaking phase transition for $c<2/3$, with $w_{12}$ and $w_{21}$ gradually deviating from each other at the critical temperature $T_{2}^{\textrm{rb}}=(2-5c^{2}/2)/(1-c^{2}/4)^{2}$ [31]. When input correlation $c$ exceeds $2/3$, a local minimum of the free energy with $w_{12}\neq w_{21}$ first emerges in the reciprocal ($w_{12}=w_{21}$) phase and it then becomes the global minimum as $T$ decreases to certain critical value strictly higher than $T_{2}^{\textrm{rb}}$. A discontinuous transition then occurs, from the reciprocal phase to a reciprocity-broken phase in which unit $2$ inhibits unit $1$ and unit $1$ excites unit $2$.

.3 Reciprocity-breaking and ideal-gas law

For multi-unit systems containing $N\geq 3$ units, our numerical optimization results and local stability analysis reveal that the optimal internal model $\bm{W}$ is reciprocal ($w_{ij}=w_{ji}$) and permutation-symmetric (all $w_{ij}$ being equal) at high temperatures, but these properties break down spontaneously as $T$ drops below the critical value [31]

[TABLE]

However, Eq. (35) does not apply for $N=2$, indicating the emergence of new collective properties for $N\geq 3$ at this continuous reciprocity-breaking phase transition. Indeed we find that the energy of the reciprocity-broken ($w_{ij}\neq w_{ji}$) optimal system exactly obeys an ideal-gas law [31]

[TABLE]

in an extended temperature interval $(T_{N}^{\textrm{iglb}},\,T_{N}^{\textrm{rb}}]$, with $T_{N}^{\textrm{iglb}}$ being the lower-bound temperature at which this ideal-gas law is deviated [Fig. 2].

One type of interaction patterns realizing the ideal-gas law is the rotation-symmetric solution with weight parameters $w_{i}$ (we require $w_{1}\leq\ldots\leq w_{N-1}$ to remove trivial degeneracies):

[TABLE]

There is cyclic dominance (CD) among the $N$ units such that unit $i+1$ is most strongly inhibited by unit $i$ and it most strongly inhibits unit $i+2$ and so on. Rotation-symmetric solutions are impossible for $N=2$, and there is only one such solution for $N=3,4$ but infinitely many degenerate solutions for $N\geq 5$ [31].

We can define a quantity $\mathcal{O}_{CD}$, which compares the most dominant global cycle of directed interactions and the reverse directed cycle of interactions, to measure the degree of cyclic dominance:

[TABLE]

where $(p_{1},p_{2},\ldots,p_{N})$ denotes a permutation of the $N$ units. For example, the weight matrix shown in the right panel of Fig. 2 has a very high value of $\mathcal{O}_{CD}=0.996$. We find through numerical computations that the ideal-gas law can be achieved by cyclic-dominant weight matrices with or without rotational symmetry and also by two-component matrices composed of excitatory and inhibitory units [31].

What is the significance of the ideal-gas law (6)? In terms of the eigenvalues $\epsilon_{i}$ of the symmetric matrix $(\bm{I}+\bm{W})^{-1}C(\bm{I}+\bm{W}^{\top})^{-1}$ we have $E=\sum_{i}\epsilon_{i}$ and $2S=\sum_{i}\ln\epsilon_{i}-\ln\textrm{Det}(\bm{C})$. The entropy is then upper-bounded by $S\leq(N/2)\ln(E/N)-(1/2)\ln\textrm{Det}(\bm{C})$ at each energy $E$, and equality holds only if all the eigenvalues $\epsilon_{i}$ are identical [31]. If this entropy upper-bound can be saturated within a continuous interval of $E$, then Eq. (6) will be achieved as a consequence of energy–information competition, and it means the most efficient utilization of energy to reach the information robustness upper-bound. An appealing statistical property is that the internal states are composed of independent components of equal magnitudes [19, 32], with $\langle x_{i}^{2}\rangle=E/N$ to achieve energy equipartition and $\langle x_{i}x_{j}\rangle=0$ for $j\neq i$ to eliminate internal pair correlations (here $\langle\cdot\rangle$ denotes averaging over all the samples). We notice a probably closely related biological fact is that the neuron pair correlations of the visual cortex are very weak even for highly correlated sensory inputs [33].

As the temperature decreases to the lower critical value $T_{N}^{\textrm{iglb}}$, it will no longer be possible for an optimal weight matrix to exactly obey Eq. (6) and at the same time keep all its complex eigenvalues above the bottomline. Then the ideal-gas law will break down and a kink of the energy function $E(T)$ will be observed. For $N=5$ and $c=0.8$, $T_{5}^{\textrm{iglb}}=0.3360$ at which the energy slope changes from $2.5$ to $2.0$ [Fig. 2]. This continuous transition is not associated with any symmetry change.

.4 Discontinuous phase transitions

The energy discontinuity at $T=0.1383$ in the example curve of Fig. 2 signifies a discontinuous phase transition. Indeed the free energy at the vicinity of this temperature has two minima and they organize into two branches with distinct energies [Fig. 2]. One branch corresponds to the one-component CD network which is almost fully inhibitory. The other branch corresponds to a two-component EI-balanced network: there is permutation symmetry within the two excitatory units $1$ and $2$ of group $g_{E}$, while the three units $3,4,5$ of group $g_{I}$ have cyclically dominant interactions within themselves and they strongly inhibit group $g_{E}$. Such a pattern is in agreement with the vital biological fact of excitation–inhibition competition and balance [34, 35, 36].

One simple measure of EI-balance is to compute the net input weights of all the individual units,

[TABLE]

A value of $\mathcal{O}_{EI}$ significantly above zero indicates a highly EI-balanced network in whcih the excitatory inputs are largely cancelled by the inhibitory inputs for most of the units. For example, $\mathcal{O}_{EI}=0.626$ for the weight matrix shown in the middle panel of Fig. 2.

After the transition to the optimal EI-balanced network, the ideal-gas law is then approximately (though not exactly) recovered at lower temperatures [Fig. 2], indicating that EI-balance is beneficial for the efficiency of energy usage. This can not be achieved if all the synaptic weights are restricted to be positive (inhibitory) [31].

The phase diagram for the system containing $N=5$ units is shown in Fig. 3. The $\beta_{1}$ and $\beta_{2}$ regions in this diagram correspond to one-component cyclic-dominant optimal networks (a subregion of $\beta_{1}$ also contains EI-balanced ones). The $\gamma_{1}$–$\gamma_{4}$ regions all correspond to different types of EI-balanced optimal networks. The phase boundaries are precisely determined by the singularities of $E(T)$, and we can also define some quantitative measures such as $\mathcal{O}_{CD}$ and $\mathcal{O}_{EI}$ to describe these different phases [31]. The ideal-gas law is exact in regions $\beta_{1}$ and $\gamma_{1}$ but it only approximately holds in $\gamma_{2}$–$\gamma_{4}$. The discontinuous transition from $\gamma_{2}$ to $\gamma_{3}$ is associated with breaking the internal rotational symmetry of group $g_{I}$ and it leads to functional differentiation of the three units $3,4,5$ (see Ref. [31] for more details). The phase diagram for the smallest multi-unit system ($N=3$) is qualitatively similar to Fig. 3; and we also report some numerical results on the larger system of $N=10$ units to demonstrate the generality of our conclusions [31].

.5 Discussion

Our theoretical work attempted to consider energy and function simultaneously in a single model, and revealed that the trade-off between energy reduction and information robustness can induce spontaneous breaking of permutation symmetry in the lateral predictive coding system and drive the formation of cyclic dominance and the sudden emergence of excitation–inhibition balance among the network units. We discovered an ideal-gas law (6) of internally representing the sensory inputs as equal-magnitude independent components. Our theoretical results are in qualitative agreement with some important biological facts such as non-reciprocal interactions, excitation–inhibition balance and weak internal pair correlations. They may also have algorithmic implications.

For simplicity we have assumed the input correlation matrix $\bm{C}$ is uniform with a single parameter $c$. A natural extension of the present work is to study the ideal-gas law (6) for heterogeneous sensory inputs. Our ongoing analytical and numerical computations confirm that Eq. (6) and continuous and discontinuous phase transition will also be observed for non-uniform matrices $\bm{C}$ [31, 39].

The existence of discontinuous phase transitions and the associated multiple free energy minima may render the optimal weight matrices hard to acquire through local and gradual Hebbian learning processes [12, 37, 38]. This issue waits to be further explored. It is also very interesting to study phase transitions in nonlinear LPC models with energy cost being a sum $\sum_{i}|x_{i}|$ of absolute values and internal prediction being rectified linear {\rm max}\bigl{(}0,\sum_{i\neq j}w_{ij}x_{j}\bigr{)} or another more complicated nonlinear function [40]. We expect to observe rich continuous and discontinuous phase transitions in these more realistic models too. These nonlinear LPC systems may lead to enhanced sensitivity to fine details in typical input vectors.

Acknowledgements.

The following funding supports are acknowledged: National Natural Science Foundation of China Grants No. 11421063, No. 11747601 and No. 12247104; National Innovation Institute of Defense Technology Grant No. 22TQ0904ZT01025. Numerical simulations were carried out at the HPC cluster of ITP-CAS and also at the BSCC-A3 platform of the National Supercomputer Center in Beijing with the help of the TRNG random number generators [41].

.6 On the definition of entropy

Suppose there is noise in the mapping from sensory input $\bm{s}=(s_{1},\ldots,s_{N})^{\top}$ to the internal representation $\bm{x}=(x_{1},\ldots,x_{N})^{\top}$,

[TABLE]

where $\xi_{i}$ for $i=1,\ldots,N$ are independent Gaussian random variables with mean zero and variance $\sigma^{2}$. For the linear lateral predictive coding (LPC) system, the deterministic function $g_{i}(\bm{s})$ is

[TABLE]

where $\bm{I}$ is the identity matrix with $I_{ii}=1$ and $I_{ij}=0$ for any $i\neq j$. The mutual information between $\bm{x}$ and $\bm{s}$ is

[TABLE]

where $H[\bm{x}]$ and $H[\bm{s}]$ are the entropies of internal states $\bm{x}$ and sensory inputs $\bm{s}$, respectively, and $H[\,\bm{x}|\bm{s}\,]$ is the conditional entropy of $\bm{x}$ given $\bm{s}$. Notice that the entropy $H[\bm{s}]$ only depends on the probability distribution $P_{\textrm{input}}(\bm{s})$ of the sensory inputs $\bm{s}$, so it could be considered as a constant. The conditional entropy $H[\,\bm{x}|\bm{s}\,]$ only depends on the property of the noise vector $(\xi_{1},\ldots,\xi_{N})^{\top}$ and so it could also be regarded as a constant. Then the mutual information $\mathcal{I}[\,\bm{x};\bm{s}\,]=H[\bm{x}]+\mathcal{I}_{0}$ with $\mathcal{I}_{0}$ being a constant. To increase this mutual information we need to increase the entropy $H[\bm{x}]$ of internal states $\bm{x}$.

The marginal distribution $P_{\textrm{output}}(\bm{x})$ of the internal state $\bm{x}$ is

[TABLE]

Equation (14) is valid for the situation of small noises ($\sigma^{2}$ being sufficiently small). The entropy of $\bm{x}$ is then

[TABLE]

Again because $H[\bm{s}]$ does not depend on the weight matrix $\bm{W}$, we see that H[\bm{x}]=-\ln\bigl{|}\textrm{Det}(\bm{I}+\bm{W})\bigr{|}+H_{0} with $H_{0}$ being a constant. For our LPC problem, the real parts of all the eigenvalues of $\bm{I}+\bm{W}$ are required to be positive. The determinant of $\bm{I}+\bm{W}$ then must be positive. This explains Eq. (3) of the main text.

.7 The two-unit system ($N=2$)

When there are only two units, the steady-state relationship between the sensory input $\bm{s}$ and the internal state $\bm{x}$ is

[TABLE]

The mean energy of the system is

[TABLE]

Here $\langle\cdot\rangle$ denotes averaging over all the input samples, and $c_{ij}$ is the element of the input correlation matrix $\bm{C}$:

[TABLE]

The entropy of the system is

[TABLE]

which only depends on the product of the two synaptic weights.

We can determine the minimum value of energy $E$ at each fixed value of $w_{12}\,w_{21}$. Some example curves of $E$ versus $w_{12}\,w_{21}$ are shown in Fig. 4 for statistically symmetric inputs with $c_{11}=c_{12}=1$, and the corresponding synaptic weights $w_{12}$ and $w_{21}$ are shown in Fig. 4. We see that $w_{12}\neq w_{21}$ when the value $w_{12}\,w_{21}$ is below certain threshold value, at each fixed value $c_{12}$ of the input correlation.

The minimum value of the free energy $F$ is determined at each value of the temperature $T$ by choosing the optimal value of $w_{12}\,w_{21}$. Then the energy $E$ as a function of $T$ can be obtained. Some example results are shown in Fig. 5 for symmetric inputs ($c_{11}=c_{22}=1$).

If the synaptic weights are restricted to be inhibitory ($w_{ij}\geq 0$), then depending on the input pair correlation $c_{12}$, there are three different outcomes [Fig. 5]. If $c_{12}$ is close to unity (say $0.95$), synaptic reciprocity $w_{12}=w_{21}$ is always kept and both $E$ and $S$ are smooth functions of $T$. At lower values of $c_{12}$ (e.g., $0.80$), interaction reciprocity breaks down discontinuously at certain critical value of $T$, with an associated sudden drop in the energy $E$ and in the entropy $S$. At still lower values of $c$ (e.g., $0.60$), the breaking of reciprocity becomes a continuous transition, and both $E$ and $S$ are continuous at the transition point but their first-order derivatives are discontinuous. The minimum-energy network will be reached when one of the synaptic weights vanishes. The minimum energy can not be arbitrarily small, because the input to one unit (say $s_{1}$) is not cancelled by the output $x_{2}$ of the other unit.

If the synaptic weights are allowed to take negative (excitatory) values, the mean energy of the network at low temperature values can be remarkably reduced [Fig. 5]. The optimal strategy to achieve sufficiently low energy $E$ appears to be that one unit inhibits and is excited by the other unit. There is either a continuous or a discontinuous phase transition for $c\leq\frac{2}{3}$ and $c>\frac{2}{3}$, respectively. The phase diagram for this case of arbitrary weights is reported in Figure 1 of the main text.

When $c_{11}\neq c_{22}$ there is no permutation symmetry to break and the continuous phase transition is absent. We find that the discontinuous phase transition is still possible for this non-symmetric case. Setting $c_{22}=1$, Fig. 6 shows the region of the $c_{11}$–$c_{12}$ plane in which the discontinuous phase transition can be observed at certain critical temperature.

.8 Local stability of the reciprocal and permutation-symmetric solution

For symmetric input vectors $\bm{s}$ with $c_{ii}\equiv\langle s_{i}^{2}\rangle=1$ and $c_{ij}\equiv\langle s_{i}s_{j}\rangle=c$ for any $i\neq j$, the corresponding symmetric internal model is the permutation-symmetric matrix $\bm{W}$ with all its non-diagonal elements being equal to $w$. Reciprocity obviously holds in this internal model ($w_{ij}=w_{ji}$ for all pairs). We now determine the value of this uniform weight $w$ and study the local stability of this reciprocal and permutation-symmetric solution.

First, notice that the homogeneous correlation matrix $\bm{C}$ could be written as

[TABLE]

from which it is easy to see that $\bm{C}$ has one eigenvalue $1+(N-1)c$ with the corresponding eigenvector \bm{u}_{0}=\frac{1}{\sqrt{N}}\bigl{(}1,1,\ldots,1\bigr{)}^{\top} and $N-1$ degenerate eigenvalues $(1-c)$ and the corresponding $N-1$ eigenvectors which are orthogonal to $\bm{u}_{0}$ and to each other. The symmetric matrix $\bm{I}+\bm{W}$ can be diagonalized in the same way. Since all the eigenvalues of the correlation matrix must be non-negative, so $c\in[-\frac{1}{N-1},\;1]$. To ensure that the short-time predictive coding dynamics be convergent, the eigenvalues of $\bm{I}+\bm{W}$ are required to be positive, which means $w\in(-\frac{1}{N-1},\;1)$. When the matrix $\bm{W}$ is permutation-symmetric, the entropy $S$ is expressed as

[TABLE]

Let us define an auxiliary symmetric matrix $\bm{M}$ for a general weight matrix $\bm{W}$ as

[TABLE]

In the special case of $\bm{W}$ being permutation-symmetric, it is easy to check that

[TABLE]

and then we obtain the following expression for the energy:

[TABLE]

By minimizing the free energy $F=E-TS$ with respect to the weight parameter $w$, we obtain the following relationship between $w$ and the temperature $T$:

[TABLE]

The local stability of this symmetric solution is studied by assuming $w_{ij}=w+\varepsilon_{ij}$ with $\varepsilon_{ij}$ being a tiny perturbation of weight $w_{ij}$. Denote by $\bm{\varepsilon}$ the perturbation matrix whose diagonal elements are all being exactly zero, we have

[TABLE]

Then the energy deviation $\Delta E$ up to the second order of $\bm{\varepsilon}$ is

[TABLE]

The first and second derivatives of the entropy are, respectively,

[TABLE]

and therefore the deviation of entropy up to second order is

[TABLE]

The first-order deviation ($\Delta F^{(1)}$) of the free energy is then

[TABLE]

If this first-order deviation is identical to zero for all perturbation matrices $\bm{\varepsilon}$ which do not violate the eigenvalue constraints of $\bm{W}$, then the matrix $\bm{W}$ must be an extremal point of the free energy $F$. The permutation-symmetric matrix with $w$ satisfying Eq. (25) is such an extremal solution.

The second-order derivation ($\Delta F^{(2)}$) of the free energy is expressed as

[TABLE]

where $i\neq j$ and $k\neq l$ and the element $L_{ij,kl}$ of the symmetric Hessian matrix $\bm{L}$ is computed according to

[TABLE]

When the minimum eigenvalue of this $N(N-1)\times N(N-1)$ symmetric matrix $\bm{L}$ changes from being positive to being negative, the synaptic matrix $\bm{W}$ will no longer be a locally stable solution of the free energy minimization problem. For the permutation-symmetric matrix $\bm{W}$ the Hessian matrix is easy to compute through Eq. (33). Our numerical results for the cases of $N=2,3,4,5,10$ are shown in Fig. 7. This reciprocity-breaking (and also permutation-symmetry-breaking) temperature is denoted as $T_{N}^{\textrm{rb}}$ in the main text.

When $N=2$ the Hessian matrix $\bm{L}$ is a $2\times 2$ matrix, and it is straightforward to obtain the analytical expression of $T_{2}^{\textrm{rb}}$ as

[TABLE]

which is in full agreement with the numerical results obtained at $N=2$.

Figure 7 demonstrates that the curve of $T_{2}^{\textrm{rb}}$ for $N=2$ is different from the corresponding curves for $N\geq 3$. After a careful examination of the numerically obtained eigenvectors of the Hessian matrices $\bm{L}$ for $N=2$ and $N\geq 3$, we realize that the systems containing $N\geq 3$ units have more ways to break the permutation symmetry and reciprocity than the two-unit system. For example, one additional symmetry for $N\geq 3$ is the rotation symmetry among the $N$ units. Our numerical results suggest that the permutation symmetry of the system starts to be unstable at $T_{N}^{\textrm{rb}}$ and it may evolve to an optimal solution with rotation symmetry among the $N$ units. This is easy to understand: a rotation-symmetric weight matrix $\bm{W}$ could be regarded as a long wavelength perturbation of a permutation-symmetric one.

Based on this insight, we compute the minimum eigenvalue of the Hessian matrix $\bm{L}$ corresponding to a rotation-symmetric perturbation vector $\{\varepsilon_{ij}\}$. We then obtain an explicit expression of the critical temperature $T_{N}^{\textrm{rb}}$ for $N\geq 3$ as

[TABLE]

by locating the point at which this eigenvalue changes sign. Another way of deriving this expression will be described in the next section. This analytical formula is in full agreement with our numerical results in Fig. 7. The $N\rightarrow\infty$ limit of this expression, $T_{N\rightarrow\infty}^{\textrm{rb}}=2(1-c)$, is also shown in this figure.

.9 On the ideal-gas law, cyclic-dominance, and excitation–inhibition balance

The ideal-gas law $E=\frac{N}{2}T$, possible only for $N\geq 3$, was first observed through numerical computations. It was found to hold true in an extended temperature range immediately after interaction reciprocity (and permutation symmetry) is broken. We now confirm these numerical results through rigorous analytical derivations.

First, as noticed in the main text, the energy and entropy can be computed as

[TABLE]

with $\epsilon_{1},\ldots,\epsilon_{N}$ being the $N$ non-negative eigenvalues of the symmetric matrix $\bm{M}$ defined in Eq. (22). By applying the following sum–product inequality

[TABLE]

we obtain the following upper bound for the entropy,

[TABLE]

This upper bound of $S$ will be attained only if all the eigenvalues $\epsilon_{i}$ take the same value $E/N$. If this is the case for a continuous range of $E$, then by applying $T=\frac{{\rm d}E}{{\rm d}S}$ as required by the minimization of the free energy $F$, we get the ideal-gas law $E=\frac{N}{2}T$.

From the ideal-gas law we can easily derive that

[TABLE]

which means that the general form of the matrix $\bm{I}+\bm{W}$ is

[TABLE]

where $\bm{U}$ is certain real orthogonal matrix satisfying $\bm{U}\bm{U}^{\top}=\bm{U}^{\top}\bm{U}=\bm{I}$, and the square-root of $\bm{C}$ is expressed as

[TABLE]

Because the diagonal elements of $\bm{W}$ is identical to zero, we see from Eq. (40) that each column of the orthogonal matrix $\bm{U}$ needs to satisfy the additonal condition of

[TABLE]

If the temperature $T$ is too large, it is impossible to satisfy these $N$ conditions (42) by any orthogonal matrix $\bm{U}$. When $T$ is lowered to the critical value $T_{N}^{\textrm{rb}}$ as given by Eq. (35), we find that these conditions are uniquely satisfied by the orthogonal matrix $\bm{U}=\bm{I}$. At this critical value, the matrix $\bm{W}$ is still permutation-symmetric and all its non-diagonal elements are

[TABLE]

For the ideal-gas law to hold at $T<T_{N}^{\textrm{rb}}$, the constraints (42) require that the $i$-th column vector of $\bm{U}$ should be staying on an $(N-1)$-dimensional ellipsoid on a $N$-dimensional hypersphere of unit radius. These $N$ constraints are compatible with the constraint of $\bm{U}$ being an orthogonal matrix, as long as $N\geq 3$. To demonstrate this most clearly, we consider the type of weight matrices with rotational symmetry.

The general form of such a matrix is

[TABLE]

where $w_{0}\equiv 1$ and $w_{1},\ldots,w_{N-1}$ are the other weight parameters. Let us define $N$ discrete frequencies as $\hat{\omega}_{j}\equiv 2\pi j/N$ for any $j=0,\ldots,N-1$ (so $\hat{\omega}_{0}=0,\;\hat{\omega}_{1}=2\pi/N,\;\ldots$). Then it is easy to check that the $j$-th eigenvalue $\lambda_{j}$ of $\bm{I}+\bm{W}$ is $\lambda_{j}=\sum_{k=0}^{N-1}w_{k}e^{ik\hat{\omega}_{j}}$ for $j=0,\ldots,N-1$, and the corresponding complex eigenvector is $\bm{r}_{j}=(1,e^{i\hat{\omega}_{j}},e^{2i\hat{\omega}_{j}},\ldots,e^{(N-1)i\hat{\omega}_{j}})^{\top}/\sqrt{N}$. It is easy to check that $\bm{r}_{j}^{\top}\bm{r}_{j}^{*}=1$ and $\bm{r}_{i}^{\top}\bm{r}_{j}^{*}=0$ for $j\neq i$. The relationship between the eigenvalues $\{\lambda_{j}\}$ and the weight parameters $\{w_{j}\}$ can also be written as

[TABLE]

If $N$ is odd so that $N=1+2N_{h}$, then one of the eigenvalues ($\lambda_{0}$) is real, and the remaining $(N-1)$ eigenvalues form $N_{h}$ complementary pairs, namely $\lambda_{j}=\lambda_{N-j}^{*}$ for $j=1,\ldots,N_{h}$. Because of the fact $e^{-i\hat{\omega}_{j}}=e^{i\hat{\omega}_{N-j}}$, we can check from equation (45) that all the weight parameters $w_{j}$ must be real-valued. If $N$ is even with $N=2N_{h}$, then two of the eigenvalues ($\lambda_{0}$ and $\lambda_{N_{h}}$) are real-valued, while the remaining $N-2$ eigenvalues form $N_{h}-1$ complementary pairs ($\lambda_{j}=\lambda_{N-j}^{*}$ for $j=1,\ldots,N_{h}-1$). Notice that $e^{-i\hat{\omega}_{N_{h}}}=-1$. We can also check that all the weight parameters as determined by Eq. (45) are again real-valued for $N$ being even.

The mean energy $E$ and the entropy $S$ for this type of matrices are, respectively,

[TABLE]

Using the sum–product inequality for the eigenvalues $|\lambda_{i}|$ (see Eq. (37)), we obtain that

[TABLE]

The equality holds only if $|\lambda_{j}|^{2}$ for all $j=1,\ldots,N-1$ are equal to each other and that

[TABLE]

for which the temperature is

[TABLE]

and therefore $E=\frac{N}{2}T$.

Let us continue to investigate the existence and the degeneracy of the rotation-symmetric matrices satisfying the ideal-gas law. There are multiple constraints on the $N$ eigenvalues $\lambda_{i}$, including Eq. (49), the conjugate property $\lambda_{j}=\lambda_{N-j}^{*}$, the condition of

[TABLE]

as inferred from $w_{0}\equiv 1$, and furthermore, that the real part of every eigenvalue $\lambda_{j}$ must be positive to ensure the convergence of the predictive coding dynamics.

First consider the limiting case of $w_{j}=w$ for all $j\geq 1$, which corresponds to $\lambda_{0}=1+(N-1)w$ and $\lambda_{j}=1-w$ for $j\geq 1$. We get from Eqs. (49) and (50) that the only solution is

[TABLE]

This point $(w,T)$ is nothing but the critical point at which the permutation-symmetric solution becomes unstable [see Eqs. (35) and (43)].

Notice that for $N=2$ it is impossible to have rotational symmetry when reciprocity is broken ($w_{12}\neq w_{21}$), so the analysis of this section does not hold. The instability of the permutation symmetry in $N=2$ is not towards rotational symmetry, but towards one unit dominating the other one.

For $N=3$, let us denote $\lambda_{1}=a_{1}+ib_{1}$ and $\lambda_{2}=a_{1}-ib_{1}$. Then we have

[TABLE]

We see that a unique rotation-symmetric solution exists for $T\in(\frac{2(1+2c)}{9},T_{c}^{\textrm{rb}})$. There is no non-trivial degeneracy within the rotation-symmetric solutions.

For $N=4$, we have two real eigenvalues $\lambda_{0}$ and $\lambda_{2}$ and two complementary complex eigenvalues $\lambda_{1,3}=a_{1}\pm ib_{1}$. The solution is determined by

[TABLE]

At each valid temperature $T$, there is only one solution of $\lambda_{0}$, and the values of $\lambda_{2}$, $a_{1}$ and $b_{1}$ are then all uniquely determined. So there is also no non-trivial degeneracy within the rotation-symmetric solutions.

For $N=5$, we have one real-valued eigenvalue $\lambda_{0}$, and two pairs of complex eigenvalues: $\lambda_{1,4}=a_{1}\pm ib_{1}$ and $\lambda_{2,3}=a_{2}\pm ib_{2}$. To eliminate the trivial symmetry we may require that $a_{1}\geq a_{2}$. Then we have

[TABLE]

Given the temperature $T$, the value of $\lambda_{0}$ is uniquely determined. The sum $a_{1}+a_{2}$ is then constrained to a fixed value but the individual values of $a_{1}$ and $a_{2}$ are not strictly constrained. Therefore there is an infinite degrees of degeneracy within the rotation-symmetric solutions. This same conclusion holds for all $N\geq 5$.

Rotation-symmetric weight matrices are special examples of more general cyclic-dominant (CD) matrices. The concept of cyclic-dominance is very easy to understand intuitively. Roughly speaking, it means one unit (say $p_{1}$) of the system is most strongly inhibited by another unit (say $p_{2}$) while $p_{1}$ only weakly affects $p_{2}$, and unit $p_{2}$ is most strongly inhibited by yet another unit $p_{3}$ while $p_{2}$ only weakly affects $p_{3}$, …, and finally unit $p_{N}$ is most strongly inhibited by unit $p_{1}$ while $p_{N}$ only weakly affects $p_{1}$. The directed cycle $p_{1}\leftarrow p_{2}\leftarrow p_{3}\leftarrow\ldots\leftarrow p_{N}\leftarrow p_{1}$ is a global closed loop of dominance. It is relatively easy to quantitatively measure the degree of cyclic-dominance in the system. For example, we may design an order parameter $\mathcal{O}_{CD}$ as

[TABLE]

where $(p_{1},p_{2},\ldots,p_{N})$ denotes a permutation of the indices of the $N$ units. We see that $\mathcal{O}_{CD}$ compares the most dominant cycle of directed interactions among $N$ units and the reverse directed cycle of interactions. If there is no or only very weak cyclic dominance then $\mathcal{O}_{CD}\approx 0$; if there is strong cyclic dominance (the limiting case being the $N$ units form a single cycle of one-way interactions but no reverse interactions), then $\mathcal{O}_{CD}$ will take a value close to unity. The order parameter $\mathcal{O}_{CD}$ for the five units shown in the right panel of Fig. 2(b) of the main text is $0.996$.

Besides rotation-symmetric solutions there are other types of matrices, such as cyclic-dominant matrices without the property of rotational symmetry and different forms of EI-balanced networks (some concrete examples are given in the next sections). This will bring in additional degrees of degeneracy within the optimal solution space, further increasing the flexibility of achieving the ideal-gas law.

The units in a cyclic-dominant optimal weight matrices $\bm{W}$ may be regarded as forming only one component in the global level. We also encounter EI-balanced optimal weight matrices in our numerical computations. Roughly speaking, an EI-balanced weight matrix $\bm{W}$ are formed by two components of units. Group $g_{I}$ contains inhibitory units $i$ whose output synaptic weights $w_{ki}$ to other units $k$ are mostly inhibitory (i.e., $w_{ki}>0$); the other group $g_{E}$ contains excitatory units $j$ whose output synaptic weights $w_{kj}$ to other units $k$ are mostly excitatory (i.e., $w_{kj}<0$). An example of EI-balanced optimal matrix is shown in Fig. 2(b) of the main text.

One quantitative measure of EI-balance is to compute the net synaptic weights on all the $N$ units,

[TABLE]

Notice that if the sum of input weights to unit $i$ is $\sum_{j\neq i}w_{ij}\approx 0$, then this unit is EI-balanced and it contributes a value $\approx 1$ to the order parameter $\mathcal{O}_{EI}$; while if all the input weights to unit $i$ are of the same sign, then it contributes a term [math] to the order parameter. The EI-balanced optimal matrix of Fig. 2(b) has a value of $\mathcal{O}_{EI}=0.626$.

Another more coarse-grained measure of EI-balance is simply to compare the number $M_{+}$ of postive synaptic weights and the number $M_{-}$ of negative synaptic weights:

[TABLE]

The EI-balanced matrix of Fig. 2(b) attains the maximum value $1.0$ in terms of this order parameter.

.10 Numerical optimization methods

Analytical solutions of the free energy minimization problem are hard to achieve for systems containing $N\geq 3$ units, especially at low temperatures. We resort to numerical optimization methods to solve this problem. Two simulated annealing search processes are implemented to construct optimal weight matrices. One process runs at a fixed temperature value $T$ and aims at reaching the global minimum of the free energy $F$. The other works at a fixed entropy value $S$ and attempts to reach the minimum energy $E$ that is compatible with this entropy.

We introduce an annealing parameter $\kappa$ and assign it a small initial value (say $\kappa=10^{2}$). Then we run $n_{0}$ (e.g., $n_{0}=10^{5}$) stochastic trials to update $\bm{W}$ at each epoch of the optimization process. After one such epoch of trials, the next epoch then follows with a slightly increased value of $\kappa$ (such as $\kappa\leftarrow 1.01\times\kappa$), initializing $\bm{W}$ by the so-far reached best weight matrix in this process. This annealing process is continued for hundreds of epochs until $\kappa$ reaches a high value (say $10^{8}$). The tentative optimal weight matrix is then further refined (if possible) by running greedy descent at $\kappa=\infty$ for $10\times n_{0}$ trials, and the final matrix $\bm{W}$ is returned as one candidate optimal solution. To make sure that the global minimum of free energy (at fixed $T$) or energy (at fixed $S$) has been achieved, we repeat this whole $\kappa$-annealing process a number of times starting from different initial conditions and random number seeds and check whether the outcomes have converged to the same minimum value.

Temperature $T$ being fixed

If the temperature $T$ is fixed during the whole simulation process, in each elementary trial we propose a new weight matrix $\bm{W}^{\prime}$ based on the current weight matrix $\bm{W}$ according to the following updating rules:

1. (a):

Single-element perturbation. The current weight matrix $\bm{W}$ is copied to the new matrix $\bm{W}^{\prime}$ but with a single randomly chosen element $w_{ij}$ being changed to $w_{ij}^{\prime}=w_{ij}+\delta_{ij}$. The deviation $\delta_{ij}$ is uniformly distributed in $(-\Delta_{1},\Delta_{1})$. 2. (b):

All-elements perturbation. The current weight matrix $\bm{W}$ is copied to the new matrix $\bm{W}^{\prime}$ but with every element $w_{ij}$ being changed to $w_{ij}^{\prime}=w_{ij}+\delta_{ij}$. The deviation $\delta_{ij}$ is an independently and uniformly distributed random number in the range $(-\Delta_{2},\Delta_{2})$.

The eigenvalues of the matrix $\bm{I}+\bm{W}^{\prime}$ are first computed. If the real part of an eigenvalue is below a bottomline value (fixed to be $10^{-5}$ in this work), then $\bm{W}^{\prime}$ is considered as invalid and the proposed update is rejected. If the eigenvalue bottomline constraint is satisfied by all the eigenvalues, then the free energy difference $\delta F$ induced by the proposed weight matrix $\bm{W}^{\prime}$ is computed. If $\delta F\leq 0$ then $\bm{W}^{\prime}$ is copied back to $\bm{W}$ (the proposed update being accepted). If $\delta F>0$, then the proposed $\bm{W}^{\prime}$ is accepted only with probability $e^{-\kappa\delta F}$ ($\bm{W}$ is replaced by $\bm{W}^{\prime}$) and it is rejected with the remaining probability $1-e^{-\kappa\delta F}$ (the old $\bm{W}$ is kept).

We choose single-element and all-elements perturbations with equal probability $1/2$ at each elementary trial. To improve the search efficiency, both $\Delta_{1}$ and $\Delta_{2}$ are dynamically adjusted so that the success rate of weight matrix change will be approximately $0.5$. We also set the lower bound $10^{-6}$ to the values of $\Delta_{1}$ and $\Delta_{2}$ to avoid the simulation process being trapped in a local minimum. If this lower bound is reach by (say) $\Delta_{1}$ in one epoch, we reset its value to $\Delta_{1}=10^{-4}$ in the next epoch to further increase the chance of escaping a local minimum of $F$.

Entropy $S$ being fixed

If the whole simulation process is run under the constraint of fixed entropy, then in each elementary trial a new weight matrix $\bm{W}^{\prime}$ is generated by modifying a single row or a single column of the old matrix $\bm{W}$. Consider a row modification which generates the $i$-th row of $\bm{W}^{\prime}$ as

[TABLE]

The constraint of fixed entropy $S$ means the following linear constraint on the $N-1$ perturbations $\delta_{ik}$:

[TABLE]

We can first generate $N-1$ independent deviations $\delta_{ik}$ by sampling in the range $(-\Delta_{1},\Delta_{1})$, and then add a correction term $\delta_{0}$ to all these deviations to satisfy the constraint (61).

Similarly, if a column $j$ of the new matrix $\bm{W}^{\prime}$ differs from the old matrix $\bm{W}$ by $w_{kj}^{\prime}=w_{kj}+\delta_{kj}$, then the deviations $\delta_{kj}$ shall satisfy the following constraint:

[TABLE]

After a new matrix $\bm{W}^{\prime}$ has been proposed by row modification (with probability one-half) or by column modification (again with probability one-half), the eigenvalues of the matrix $\bm{I}+\bm{W}^{\prime}$ are first computed. If the real part of an eigenvalue is below the bottomline value ($10^{-5}$), then the proposed update is rejected. If the eigenvalue bottomline constraint is satisfied by all the eigenvalues, then the energy difference $\delta E$ induced by the proposed weight matrix $\bm{W}^{\prime}$ is computed. If $\delta E\leq 0$ then $\bm{W}$ is replaced by $\bm{W}^{\prime}$. If $\delta E>0$, then the proposed $\bm{W}^{\prime}$ is accepted with probability $e^{-\kappa\delta E}$ and it is rejected with the remaining probability $1-e^{-\kappa\delta E}$.

.11 Some example matrices for $N=5$

The phase diagram for $N=5$ is shown in Fig. 3 of the main text. Here we offer some additional concrete samples of optimal weight matrices to help appreciate this phase diagram. The different phase regions are discussed by decreasing the temperature $T$ along a vertical line of fixed pair correlation $c$. Here we mainly consider $c=0.6$, for which the relationship between $E$ and $S$ is shown in Fig. 8.

Phase region $\alpha$

For example, at the point $c=0.6$ and $T=1.6507$ (slightly above the critical temperature $1.5305$), the permutation-symmetric matrix is

[TABLE]

with entropy $S=0.6$ and energy $E=3.9030$, and the eigenvalues of $\bm{I}+\bm{W}$ being $\{2.1602,\;0.7100\pm 8\times 10^{-7},\;0.7100\pm 2\times 10^{-7}\}$.

Phase region $\beta_{1}$

We show two degenerate optimal matrices here for the point $c=0.6$, $T=0.45152$ in phase region $\beta_{1}$:

[TABLE]

These two solutions have the same entropy $S=-2.50$ and the same energy $E=1.1288$. The left one is a one-component system with cyclic dominant (CD) interactions but there is no rotation symmetry or other obvious symmetry, and the eigenvalues of $\bm{I}+\bm{W}$ are $\{3.8092,\;0.5785\pm 1.2062i,\;0.0169\pm 1.3367i\}$. The right network is an EI-balanced two-component system with group $g_{E}$ containing units $1$ and $2$ and group $g_{I}$ containing the remaining three units, and the eigenvalues are $\{1.3427,\;1.8103\pm 1.3213i,\;0.0184\pm 1.3439i\}$. There is no permutation symmetry within $g_{E}$, and there is cyclic dominance but no rotation symmetry in $g_{I}$.

Phase range $\beta_{2}$

For $c=0.6$, the EI-balanced solution mentioned above reaches the eigenvalue bottomline at a critical temperature value $0.3872$ and its free energy starts to be higher than that of the CD solution. The CD solution will reach the eigenvalue bottomline at a much lower critical temperature $T=0.2720$, at which its energy susceptibility drops from $5/2$ to $2.0$. This transition is completely due to eigenvalue constraint; there is no any qualitative change in the symmetry of the solution. Here is an example matrix in this phase, obtained at $c=0.6$ and $T=0.2303$:

[TABLE]

whose entropy is $S=-4.1$ and energy is $E=0.56958$ (the eigenvalues of $\bm{I}+\bm{W}$ are $\{4.99996,\;10^{-5}\pm 1.8639i,\;10^{-5}\pm 1.8637i\}$. This is a CD matrix with all the complex eigenvalues touching the bottomline.

Phase region $\gamma_{2}$

At $c=0.6$, although the EI solutions are suppressed by the CD solutions in the phase space region $\beta_{2}$, this suppression will be reversed at a critical temperature $T=0.2195$, and a discontinuous transition between the CD and EI solutions occurs here. The EI solutions at the lower-temperature phase region $\gamma_{2}$ appears to be permutation-symmetric within group $g_{E}$ and cyclic-dominant with rotation symmetry within group $g_{I}$. For example, the EI solution at $T=0.2138$ (immediately after the discontinuous transition) is

[TABLE]

with entropy $S=-4.35$ and $E=0.5414$. The eigenvalues of $\bm{I}+\bm{W}$ are $\{2.0086,1.4957\pm 2.8409i,10^{-5}\pm 1.9345i\}$.

Phase region $\gamma_{3}$

As the temperature further decreases to another critical value $T=0.09505$ (with $c=0.6$), another discontinuous transition occurs to the optimal EI-balanced solution, caused by the breaking of internal rotation symmetry among the three units of group $g_{I}$. This transition is associated with a small drop of the entropy from $-6.410$ to $-6.417$ and a small drop of the mean energy from $0.2397$ to $0.2390$. To demonstrate the different symmetries of these two types of EI-balanced solutions, we show the two corresponding solutions at the same entropy value $S=-6.414$ and the same energy value $E=0.2394$:

[TABLE]

Notice that the rotation-symmetry breaking of the group $g_{I}$ does not involve all the five units of the system; it is an abrupt change only in the subnetwork $g_{I}$ of the system. Notice also that the permutation symmetry between the two units $1$ and $2$ of group $g_{E}$ is intact in the EI-balanced phase space region $\gamma_{3}$.

Phase region $\gamma_{4}$

As the temperature $T$ further decreases, the EI-balanced optimal solutions of the phase region $\gamma_{3}$ will reach the eigenvalue bottomline, leading to a kink in the mean energy. This continuous transition occurs at $T=0.0383$ for $c=0.6$ (with entropy $S\approx-8.73$). Notice that there is no any qualitative change in the symmetry of the EI-balanced matrix, only that all the four complex eigenvalues touch the bottomline.

As the temperature decreases even further, the permutation symmetry between the two units of group $g_{E}$ will also break. This breaking of internal symmetry will be associated with a very weak continuous transition and a tiny kink in the mean energy. At very low temperatures, then there are no internal symmetries within both group $g_{E}$ and $g_{I}$. To give one such example, we show the matrix at a very low temperature $T=0.02652$ with entropy $S=-9.47$ and mean energy $E=0.071037$ here:

[TABLE]

The eigenvalues of $\bm{I}+\bm{W}$ are $\{4.99996,10^{-5}\pm 7.1819i,10^{-5}\pm 7.0902i\}$. Notice that cyclic dominance still exist in group $g_{I}$ and the three units of this group are all different.

Phase region $\gamma_{1}$

Here we show an optimal solution obtained at $c=0.4$ and $T=0.0978$ in this second ideal-gas region $\gamma_{1}$:

[TABLE]

The entropy is $S=-7$ and energy is $E=0.2446$. The eigenvalues of $\bm{I}+\bm{W}$ are $\{3.5236,0.6740\pm 4.9114i,0.0642\pm 3.5581i\}$. This is an EI-balanced network with two groups of units. This is no permutation symmetry between the two units $1$ and $2$ of group $g_{E}$, but there is cyclic dominance within the three different units ($3,4,5$) of group $g_{I}$.

Phase region $\delta$

This phase region is located at $c$ close to unity and $T$ relatively small. There is a discontinuous phase transition between the $\beta_{2}$ phase and the $\delta$ phase. For example at $c=0.99$ this occurs at $T=0.0612$. The optimal weight matrix at temperature $T=0.0613$ (slightly above the critical temperature) is

[TABLE]

which has entropy $S=0.631$ and energy $E=0.3210$ (the eigenvalues of the matrix $\bm{I}+\bm{W}$ are $\{4.99996,10^{-5}\pm 0.5712i,10^{-5}\pm 0.5711i\}$. It is a one-component cyclic-dominant network. Immediately after the discontinuous transition, the optimal matrix at $T=0.0608$ is

[TABLE]

whose entropy is $S=-1.871$ and energy $E=0.1679$ (the eigenvalues of $\bm{I}+\bm{W}$ are $\{4.99996,\;10^{-5}\pm 1.0680i,\;10^{-5}\pm 1.0672i\}$). Notice that unit $1$ excites units $2,3,5$ but inhibits unit $4$, while the remaining four units are all inhibitory units.

At a lower critical temperature value $T=0.03445$ there is again a discontinuous phase transition, from an optimal network of the $\delta$ phase to an optimal network of the $\gamma_{3}$ phase. At the phase transition point, the optimal network of the $\delta$ phase is

[TABLE]

and its entropy is $S=-3.50$ and energy is $E=0.0914$ (the eigenvalues of $\bm{I}+\bm{W}$ are $\{4.99996,\;10^{-5}\pm 1.6044i,\;10^{-5}\pm 1.6041i\}$. The optimal network after this transition is

[TABLE]

whose entropy is $S=-3.6$ and energy is $E=0.08796$ (the eigenvalues of $\bm{I}+\bm{W}$ are $\{0.8492,\;2.0754\pm 6.5191i,\;10^{-5}\pm 0.9596i\}$.

.12 Ideal-gas law and phase diagram for the three-unit system ($N=3$)

The ideal-gas law for general input correlation matrix

First, it is instructive to study the $N=3$ system with symmetric input with a single parameter $c$ for the pair correlations by exhaustive search. To achieve the ideal-gas law, the following condition must be satisfied by the synaptic weight matrix

[TABLE]

We consider the following general form of $\bm{I}+\bm{W}$ as

[TABLE]

We may require $c_{1}\leq c_{2}$ to eliminate the trivial degeneracy of the synaptic weight matrix. Then we obtain that

[TABLE]

To determine the weight parameters, we express

[TABLE]

with $\theta\in[\frac{\pi}{4},\frac{5\pi}{4}]$, and $x=\sqrt{\frac{2}{T}-1}\geq 0$ (so $T\leq 2$). Then we find that

[TABLE]

where $y=\frac{2c}{T}$. Similarly, the expressions for $a_{2}$ and $b_{2}$ are

[TABLE]

Notice that the following condition must be satisfied by $a_{1},b_{1},a_{2},b_{2}$:

[TABLE]

which means that

[TABLE]

The other constraints (to ensure the weight values being real) are

[TABLE]

All the eigenvalues of $\bm{I}+\bm{W}$ should have positive real part.

At each value of $\theta$ there are up to four distinct solutions of $(a_{1},a_{2},b_{1},b_{2})$. They can be determined accurately by numerical computations. The numerical results are the following.

For the EI-balanced network, numerical results indicate that the upper-bound temperature can be obtained by setting $c_{1}=c_{2}=-\frac{x}{\sqrt{2}}$. Then $b_{1}=b_{2}$ and $a_{1}=a_{2}$. We can then obtain the upper-bound temperature. For example, the upper-bound EI-balanced temperature at $c=0.35$ is $T=0.264684$, and the corresponding network is

[TABLE]

This EI-balanced system breaks the internal permutation symmetry within group $g_{E}$ (units $2$ and $3$) to achieve the ideal-gas law. On the other hand, the lower-bound temperature for the ideal-gas law for the EI-balanced network can be computed numerically by computing the eigenvalues of $\bm{W}$. The temperature upper-bound and lower-bound are shown in Fig. 9.

The other types of network realizing the ideal-gas law can all be obtained. For $N=3$ these are cyclic-dominant networks. The temperature upper-bound and lower-bound for these ideal-gas CD networks are shown in Fig. 9. By comparing Figs. 9 and 9 we see that there is a region in the phase space in which EI-balanced and CD optimal solutions coexist.

The above analysis also works for the case of $\bm{C}$ having different values:

[TABLE]

with parameters $c_{22}$, $c_{33}$, $c_{12}$, $c_{13}$, and $c_{23}$ (the value of $c_{11}$ is fixed to be $c_{11}=1$ and we require that $1\leq c_{22}\leq c_{33}$). Then we still have Eq. (77) for the weight parameters $c_{1}$ and $c_{2}$, and the equations for $a_{1}$ and $b_{1}$ are

[TABLE]

where $x_{2}=\sqrt{\frac{2c_{22}}{T}-1}$ and $y_{1}=\frac{2c_{12}}{T}$. Similarly, the expressions for $a_{2}$ and $b_{2}$ are

[TABLE]

where $x_{3}=\sqrt{\frac{2c_{33}}{T}-1}$ and $y_{2}=\frac{2c_{13}}{T}$. The other constraint is

[TABLE]

where $y_{3}=\frac{2c_{23}}{T}$.

We now examine a particular example, with correlation matrix

[TABLE]

Examining all the possible solutions of the ideal-gas-law condition (84), we find that the ideal-gas law is achievable for temperatures $T\in[0.19402,1.98262]$. At the maximum temperature $T=1.98262$ the ideal-gas matrix is

[TABLE]

In temperature range $1.70899<T<1.98262$ there are two degenerate ideal-gas matrices. For example, at $T=1.8$ they are

[TABLE]

Notice that $\bm{W}^{(1)}$ is cyclic-dominant, and $\bm{W}^{(2)}$ contains a fully inhibitory unit, a fully excitatory unit, and a mixed unit.

In the temperature range $1.42383<T<1.70899$ there are four degenerate ideal-gas matrices. For example at $T=1.6$ they are:

[TABLE]

Matrices $\bm{W}^{(1)}$ and $\bm{W}^{(3)}$ are cyclic-dominant, and $\bm{W}^{(2)}$ and $\bm{W}^{(4)}$ contain a fully inhibitory unit, a fully excitatory unit, and a mixed unit.

In the temperature range $1.1927<T<1.42383$ there are six degenerate ideal-gas matrices. For example at $T=1.3$ they are:

[TABLE]

Matrices $\bm{W}^{(1)}$ and $\bm{W}^{(3)}$ are cyclic-dominant, and $\bm{W}^{(2)}$ and $\bm{W}^{(4)}$ contain a fully inhibitory unit, a fully excitatory unit, and a mixed unit, and matrices $\bm{W}^{(5)}$ and $\bm{W}^{(6)}$ are two excitation-inhibition balanced network with one inhibitory unit and two excitatory units.

In the temperature range $0.31636<T<1.1927$ there are eight degenerate ideal-gas matrices. For example at $T=0.7$ they are:

[TABLE]

Matrices $\bm{W}^{(1)}$ and $\bm{W}^{(3)}$ are cyclic-dominant, and the other six matrices are all formed by three functionally distinct units.

The degeneracy then decreases at lower temperatures. First the matrix $\bm{W}^{(1)}$ attaches the eigenvalue bottomline at $T=0.31636$, followed by matrix $\bm{W}^{(3)}$ at $T=0.31346$, followed by $\bm{W}^{(6)}$ at $T=0.23364$, followed by $\bm{W}^{(8)}$ at $T=0.22835$, followed by $\bm{W}^{(5)}$ at $T=0.21387$, followed by $\bm{W}^{(4)}$ at $T=0.21290$, followed by $\bm{W}^{(2)}$ at $T=0.19830$, and finally followed by $\bm{W}^{(7)}$ at $T=0.19402$. At the lowest temperature $T=0.19402$ the ideal-gas law matrix is

[TABLE]

The above detailed results on the problem instance (88) demonstrate that the abundance of ideal-gas optimal matrices. These results indicate that, when the temperature $T$ is in an appropriate range, the ideal-gas law optimal solutions might be relatively easy to obtain. We expect these aspects also extend to higher-dimensional systems with $N>3$ units.

Phase diagram for uniform input correlation matrix

For the $N=3$ system with symmetric sensory inputs the phase diagram in the plane of $c$ and $T$ is shown in Fig. 10. This phase diagram is qualitatively similar to the one for $N=5$ (Fig. 3 of the main text) and it is much more richer than the phase diagram for $N=2$ (Fig. 1 of the main text), confirming the intuitive expectation that three-body and multi-body (effective) interactions are capable of producing high degrees of complexity.

We briefly describe the properties of these different phase space regions here.

1. $\alpha$:

one-component networks with permutation symmetry. 2. $\beta_{1}$:

containing both one-component networks with cyclic dominance and two-component networks with excitation-inhibition balance. The ideal-gas law is obeyed. 3. $\beta_{2}$:

one-component networks with cyclic dominance, constrained by eigenvalue bottomline. 4. $\gamma_{1}$:

two-component EI-balanced networks without internal permutation symmetry within the group of two units. The ideal-gas law is obeyed. One example network obtained at $c=0.35$ and $T=0.2643$ (slightly below the $\gamma_{2}\rightarrow\gamma_{1}$ critical temperature value $0.2644$) is

[TABLE]

whose entropy is $S=-2.87$ and energy $E=0.3965$, and the eigenvalues of $\bm{I}+\bm{W}$ are: $\{2.2177,\;0.3912\pm 2.7928i\}$. 5. $\gamma_{2}$:

two-component EI-balanced networks with internal permutation symmetry within the group of two units. One example for $c=0.35$ and $T=0.2661$ (slightly above the $\gamma_{2}\rightarrow\gamma_{1}$ critical temperature $0.2644$) is

[TABLE]

whose entropy $S=-2.86$ and energy $E=0.3992$) and the eigenvalues of $\bm{I}+\bm{W}$ are $\{2.2106,\;0.3947\pm 2.7827i\}$. 6. $\gamma_{3}$:

two-component EI-balanced networks at high values of $c$ without permutation symmetry within the group of two units. 7. $\gamma_{4}$:

two-component EI-balanced networks at low values of $c$ constrained by the eigenvalue bottomline. One example obtained at $c=0.35$ and $T=0.1404$ is

[TABLE]

whose entropy $S=-3.8$ and energy $E=0.21335$, and the eigenvalues of $\bm{I}+\bm{W}$ are $\{2.99998,\;10^{-5}\pm 3.8601i\}$. There is no permutation symmetry within the group of two units. 8. $\gamma_{5}$:

two-component EI-balanced networks at high values of $c$ constrained by the eigenvalue bottomline. There is no permutation symmetry within the group of two units. 9. $\gamma_{6}$:

two-component EI-balanced networks constrained by the eigenvalue bottomline, with permutation symmetry within the group of two units. The transition between $\gamma_{4}$ and $\gamma_{6}$ and that between $\gamma_{5}$ and $\gamma_{6}$ are both continuous. An example obtained at $c=0.35$ and $T=0.0464$ is

[TABLE]

whose entropy $S=-4.9$ and $E=0.1198$ and the eigenvalues of $\bm{I}+\bm{W}$ are $\{2.99998,\;10^{-5}\pm 6.6906i\}$. The permutation within the two excitatory units $2$ and $3$ is preserved in phase region $\gamma_{6}$. 10. $\delta$:

An intermediate phase region between $\beta_{2}$ and $\gamma_{3}$ at high values of $c$.

Our numerical computations suggest that the intermediate phase region $\delta$ only exists for $c>0.9678$. To demonstrate the discontinuous transitions between $\beta_{2}$ and $\delta$ and between $\delta$ and $\gamma_{3}$, we plot the free energy curves of these three phases for $c=0.99$ at low temperatures $T$ close to $0.09$ in Fig. 11.

There is a discontinuous transition between the $\beta_{2}$ phase and the $\delta$ phase at $T=0.09905$, at which the weight matrix changes from a purely inhibitory one with $S=0.5$ and $E=0.4300$ (the eigenvalues of $\bm{I}+\bm{W}$ are $\{2.99998,\;10^{-5}\pm 0.4496i\}$) to one with a single strong excitatory weight and five inhibitory weights whose $S=-2.1$ and $E=0.1725$ (the eigenvalues of $\bm{I}+\bm{W}$ are $\{2.99998,\;10^{-5}\pm 1.6499i\}$:

[TABLE]

This phase transition brings a large drop in the energy and in the entropy. At a lower temperature $T=0.08369$, there is another discontinuous transition, from a matrix of the $\delta$ phase with $S=-2.45$ and $E=0.1407$ (the eigenvalues of $\bm{I}+\bm{W}$ are $\{2.99998,\;10^{-5}\pm 1.9645i\}$) to an EI-balanced matrix of the $\gamma_{3}$ phase with $S=-2.6$ and $E=0.1281$ (the eigen values of $\bm{I}+\bm{W}$ are $\{0.6148,\;1.1926\pm 4.5251i\}$):

[TABLE]

This second discontinuous transition only causes a relatively small drop in the energy and in the entropy, but the change in the optimal matrix $\bm{W}$ is quite comprehensive.

.13 Some numerical results for $N=10$ and $c=0.8$

Besides $N=3$ and $N=5$, we also study systems with other number of units up to $N=10$. The phase diagram for the system with $N=4$ units has also been obtained by numerical computations. It is qualitatively similar to the phase diagrams of $N=3$ and $N=5$. Here we show some numerical results obtained on the large system of $N=10$, which reveal the competition among multiple local minima of the free energy. The input pair correlation is fixed to be $c=0.8$ for simplicity.

We confirm that there is a continuous permutation-symmetry-breaking transition at a high critical temperature $T_{c}^{(1)}=0.9489$ for $N=10$ and $c=0.8$, at which cyclic dominance starts to emerge, and the energy versus temperature relationship starts to follow the ideal-gas law $E=5T$. The optimal solutions are not unique but are degenerate. For example, we observe two distinct solutions $\bm{W}$ with the same entropy $S=-4.6$ and the same mean energy $E=1.1555$ at temperature $T=0.2311$. One of them is one-component with mostly positive synaptic weights, while the other is a two-component EI-balanced solution containing a single purely excitatory unit and nine inhibitory units.

As temperature decreases to another critical value $T_{c}^{(2)}=0.1763$, the ideal-gas law $E=5T$ finally breaks as all the complex eigenvalues of the optimal solutions reach the bottomline. The energy susceptibility then drops from $5$ to $4.5$. At temperatures slightly below this critical value the minimum free energy solution is the two-component EI-balanced one with group $g_{E}$ containing a single excitatory unit and group $g_{I}$ containing the remaining nine inhibitory units (we call this the $1e$–$9i$ type structure).

As temperature decreases to much below $T_{c}^{(2)}$, we find that two other types of EI-balanced networks, the $3e$–$7i$ type with $g_{E}$ containing three excitatory units and $g_{I}$ containing seven inhibitory units and the $5e$–$5i$ type with $g_{E}$ containing five excitatory units and $g_{I}$ containing five inhibitory units, start to compete with this $1e$–$9i$ type networks. To demonstrate this competition, we draw in Fig. 12 the free energy differences of these three types of EI-balanced solutions at the vicinity of $T=0.107$.

The free energy $f_{5}$ of the $5e$–$5i$ solution starts to be lower than the value $f_{3}$ of the $3e$–$7i$ solution as temperature is decreased below the threshold value $0.1084$, so the $5e$–$5i$ structure becomes more favorable than the $3e$–$7i$ structure. At the lower critical value $0.1077$ the $5e$–$5i$ structure also starts to be dominating over the $1e$–$9i$ structure. The $3e$–$7i$ structure also has lower free energy than the $1e$–$9i$ structure at $T<0.1070$, but it can not compete with the $5e$–$5i$ structure. We therefore conclude that there is a discontinuous phase transition at $T_{c}^{(3)}=0.1077$ between the $1e$–$9i$ and $5e$–$5i$ EI-balanced network structures. This transition is associated with a discontinuity in the entropy and in the mean energy [Fig. 12].

To have a concrete feeling of the three types of EI-balanced network structures, we list three example matrices here. First is the $1e$–$9i$ network with entropy $S=-8.35$:

[TABLE]

whose energy is $E=0.55036$ and the eigenvalues of $\bm{I}+\bm{W}$ are $\{7.083736231,\;2.916183746,\;10^{-5}\pm 1.9483i,\;10^{-5}\pm 1.9477i,\;10^{-5}\pm 1.9472i,\;10^{-5}\pm 1.9464i\}$. The $3e$-$7i$ network with entropy $S=-8.35$ is:

[TABLE]

whose energy is $E=0.55031$ and the corresponding eigenvalues are $\{2.0196,\;2.0196,\;2.9803\pm 3.4025i,\;10^{-5}\pm 1.9277i,\;10^{-5}\pm 1.9268i,\;10^{-5}\pm 1.9265i\}$. There is permutation symmetry within the three units of group $g_{E}$. The $5e$–$5i$ network with $S=-8.36$ is

[TABLE]

whose energy is $E=0.55032$ and the eigenvalues of $\bm{I}+\bm{W}$ are $\{1.9942\pm 2\times 10^{-6}i,\;1.9942\pm 2\times 10^{-5},\;1.0115\pm 4.3843i,\;10^{-5}\pm 1.9117i,\;10^{-5}\pm 1.9108i\}$. There is permutation symmetry within the five units of group $g_{E}$ and cyclic dominance within the five units of group $g_{I}$.

.14 The effect of restricting all the synaptic weights to be non-negative

Some neurons in the biological brain are purely inhibitory, and the networks formed by them will have no internal excitations. If the predictive coding system is composed purely of inhibitory units, then all the synaptic weights will be non-negative. Figure 5 has clearly demonstrated the significant effect at low temperatures of restricting the synaptic weights to be non-negative. Here we present some additional results obtained on multi-unit systems ($N=5$). To compare with Fig. 2 of the main text, we fix the input correlation to be $c=0.8$, see Fig. 13.

If the temperature is high enough the restriction of all weights being non-negative has no effect to the system’s property. The weight matrix is completely symmetric at high temperatures ($w_{ij}=w_{ji}=w>0$). A continuous transition occurs at $T=1.179$ with reciprocity breaking ($w_{ij}\neq w_{ji}$) and with the formation of cyclic dominance among any set of three units and the emergence of one global strongest inhibition cycle involving all the units. The energy susceptibility, $\chi\equiv{\rm d}E/{\rm d}T$, has a discontinuity at this transition point, and the ideal-gas law follows after this permutation-symmetry breaking transition.

As $T$ decreases to $0.336$ the susceptibility drops from $2.5$ to $2.0$ but the energy $E$ itself is continuous, indicating another continuous transition. Another signature of this transition is that the real parts of all the complex eigenvalues of $\bm{I}+\bm{W}$ touch the bottomline value $10^{-5}$. At this transition point all the synaptic weights are still positive, so the non-negative constraint does not affect this transition.

As $T$ further decreases to $0.15$ some of the weights reach zero and the mean energy $E$ starts to level off. There is then a considerable drop in $E$ at $T=0.128$ which is caused by a sudden relatively huge change in the weight matrix (see Fig. 13, $T=0.127$ and $0.130$).

As $T$ keeps decreasing, the weight matrix experiences several additional abrupt changes before reaching the completely one-directional form at very low $T$ values. The discontinuous nature of weights change is confirmed by the strong hysteresis behavior of energy at low temperatures [Fig. 13, inset]. The nearly minimum-energy weight matrix at $T=0.013$ has a clear hierarchical structure: neuron $1$ is inhibited by neuron $5$ and inhibits neuron group $\{2,3,4\}$ which in turn inhibits neuron $5$, and there are internal cyclic inhibitions within group $\{2,3,4\}$ and within group $\{1,2,5\}$. Because excitation–inhibition balance can not be realized in the system, at low temperatures the ideal-gas law has no chance of being approximately restored. The mean energy saturates to a value slightly below $0.35$ as $T$ decreases to zero.