For the manuscript titled “Self-orthogonalizing attractor neural networks emerging from the free energy principle” by Tamas Spisak and Karl Friston¶

Library imports for inline code.

Appendix 1¶

Continuous Bernoulli distribution

$x \sim \mathcal{CB}(L) \iff P(x) \propto e^{Lx}, \quad x \in [-1, 1]\subset \mathbb{R}$

Below we provide a Python implementation of the continuous Bernoulli distributions: $\mathcal{CB}(L)$, parametrized by the log odds L and adjusted to the [-1,1] interval.

Appendix 2¶

Continuous Bernoulli distribution full derivation

Derivation of the exponential form of the continuous Bernoulli distribution, parametrized by the log odds L and adjusted to the [-1,1] interval.

$P(x; L) = \frac{e^{Lx}}{\int_{-1}^1 e^{Lx}} = \frac{e^{Lx}}{ \frac{e^L - e^{-L}}{L}} = L \frac{e^{Lx}}{2sinh(L)}$

Appendix 3¶

Visualization of the likelihood function with various parameters, in python.

Appendix 4¶

Expected value of the $\mathcal{CB}$

Appendix 5¶

Conservative dynamics

Let $x ∈ ℝ^n$ denote the system’s states (internal, blanket, and external states). The drift can be decomposed into a gradient part (from a potential U) and a solenoidal part R:

$\dot{x} = −∇U(x) + R(x)$

where $R = −Rᵀ$ is antisymmetric in state‐space.

The probability density $p(x, t)$ over states evolves according to the Fokker–Planck equation:

$∂ₜp(x,t) = −∇ · [ (−∇U + R) p(x,t) ] + diffusion terms.$

To determine the stationary distribution $p_s(x)$, we set $∂ₜp(x,t) = 0$. The Fokker-Planck equation then implies that the net probability current $J_s(x) = (−∇U(x) + R(x)) p_s(x) - D∇p_s(x)$ (assuming isotropic diffusion $D$) must be divergence-free: $∇ \cdot J_s(x) = 0$. If we propose $p_s(x) = Z^{-1} \exp(−U(x)/D)$ (where $Z$ is a normalization constant), then the diffusive current $-D∇p_s(x)$ becomes $p_s(x)∇U(x)$.

Substituting this into $J_s(x)$, we get:

$J_s(x) = (−∇U(x) + R(x)) p_s(x) + p_s(x)∇U(x) = R(x) p_s(x)$.

Thus, $p_s(x) \propto \exp(−U(x)/D)$ is the stationary distribution if and only if the solenoidal component of the probability current, $J_R(x) = R(x) p_s(x)$, is itself divergence-free: $∇ \cdot (R(x) p_s(x)) = 0$.

It is a key result from the study of non-equilibrium systems under the Free Energy Principle, particularly those possessing a Markov blanket and involving conservative particles, that this condition $∇ \cdot (R(x) p_s(x)) = 0$ holds (see Friston et al., 2023; but also related: Ao (2004); Xing (2010)). The conditional independence structure imposed by the Markov blanket, along with other FEP-specific assumptions (like conservative particles), constrains the system’s dynamics such that solenoidal forces $R(x)$ do not alter the Boltzmann form of the stationary distribution. Instead, they drive persistent, divergence-free probability currents $J_R(x)$ along iso-potential surfaces.

In steady state $(∂ₜp = 0)$, only the gradient term $−∇U$ influences $p(x)$. The stationary (nonequilibrium) distribution is

$p(x) ∝ exp{ −U(x) }$,

determined solely by the symmetric (potential) part −∇U. The antisymmetric $R$ does not reweight $p(x)$; it just circulates probability around isocontours of $U$.

In most NESS systems, antisymmetric flows do alter the stationary measure. Under particular‐partition (Markov‐blanket) constraints, however, the internal–external factorization ensures that solenoidal (antisymmetric) flows remain divergence‐free, leaving the Boltzmann‐like steady distribution intact. This is why, in these Markov‐blanketed systems, one can have persistent solenoidal currents (nonequilibrium flows) yet preserve a stationary distribution that depends only on the symmetric part of the couplings.

Appendix 6¶

Detailed inference derivation: accuracy–complexity decomposition and $\dfrac{\partial F}{\partial b_q}$

The main text derives the inference rule (eq. (15)) via a compact local VFE argument. Here we provide the equivalent — but more detailed — derivation via the accuracy–complexity decomposition of VFE, confirming the same result.

1. Subsitute our parametrization into F

Let’s start with substituting our parametrization into eq. .

1a. Accuracy term

From the RBM marginalization (eq. ):

$E(\bm{\sigma}) = \underbrace{-b_i\sigma_i - \sum_{j\neq i}J_{ij}\sigma_i\sigma_j}_{\text{Terms with } \sigma_i} \underbrace{-\sum_{j\neq i}b_j\sigma_j - \frac{1}{2}\sum_{j,k\neq i}J_{jk}\sigma_j\sigma_k}_{\text{Terms without } \sigma_i}$

Here, $-b_i\sigma_i$ becomes constant, since $\sigma_i$ is fixed. So we get:

$P(\sigma_{\backslash i}|\sigma_i) \propto \exp\left(\sum_{j\neq i}(b_j + J_{ij}\sigma_i)\sigma_j + \frac{1}{2}\sum_{j,k\neq i}J_{jk}\sigma_j\sigma_k\right)$

Taking expectation of $\ln P(\sigma_{\backslash i}|\sigma_i)$ under $q(σ_i)$:

$\mathbb{E}q[\ln P(σ_{\backslash i}|σ_i)] = \text{const} + \sum_{j\neq i}b_jσ_j + S(b_q)\sum_{j\neq i}J_{ij}σ_j + \frac{1}{2}\sum_{j,k\neq i}J_{jk}σ_jσ_k$

Where $S(b_q) = \mathbb{E}_q[σ_i] = \coth b_q - 1/b_q$ is the expected value of the $\mathcal{CB}$, a sigmoid function of the bias (#supplementary-information-4)).

1b. Complexity term

The complexity term in eq. is simply the KL-divergence term between two $\mathcal{CB}$ distributions. For $\mathcal{CB}$ distributions:

$q(x) = \frac{b_q}{2\sinh b_q}e^{b_q x},\quad p(x) = \frac{b}{2\sinh b}e^{b x}$

• KL divergence definition:

$D_{KL} = \int_{-1}^1 q(x) \ln\frac{q(x)}{p(x)} dx = \mathbb{E}_q[\ln q(x) - \ln p(x)]$

• Expand log terms:

$\ln q(x) = \ln b_q - \ln(2\sinh b_q) + b_q x$

$\ln p(x) = \ln b - \ln(2\sinh b) + b x$

• Subtract log terms:

$\ln\frac{q(x)}{p(x)} = \ln\frac{b_q}{b} + \ln\frac{\sinh b}{\sinh b_q} + (b_q - b)x$

• Take expectation under q(x):

$D_{KL} = \ln\frac{b_q}{b} + \ln\frac{\sinh b}{\sinh b_q} + (b_q - b)\mathbb{E}_q[x]$

• Compute expectation $\mathbb{E}_q[x]$:

$\mathbb{E}_q[x] = \int_{-1}^1 x \frac{b_q e^{b_q x}}{2\sinh b_q} dx = \frac{1}{2\sinh b_q}\left[\frac{e^{b_q x}}{b_q^2}(b_q x - 1)\right]_{-1}^1$

• Evaluate at bounds:

$= \frac{1}{2\sinh b_q}\left(\frac{e^{b_q}(b_q - 1) - e^{-b_q}(-b_q - 1)}{b_q^2}\right)$

• Simplify using hyperbolic identities:

$= \frac{(b_q \cosh b_q - \sinh b_q)}{b_q^2 \sinh b_q} = \coth b_q - \frac{1}{b_q}$

• Final substitution for the complexity term:

$D_{KL} = \ln\frac{b_q \sinh b}{b \sinh b_q} + (b_q - b)\left(\coth b_q - \frac{1}{b_q}\right)$

1c. Combining the two terms

Combining the two terms, we get the following expression for the free energy:

$F = \ln\left(\frac{b_q}{b}\right) + \ln\left(\frac{\sinh(b)}{\sinh(b_q)}\right) + (b_q - b) S(b_q) - \sum_{j \ne i} \left( b_j + S(b_q) J_{ij} \right) \sigma_j - \dfrac{1}{2} \sum_{j \ne i} \sum_{k \ne i} J_{jk} \sigma_j \sigma_k + C$

where C denotes all constants in the equation that are independent of $\sigma or b_q$.

2. Free Energy partial derivative calculation

• First, we differentiate the log terms:

$\frac{\partial}{\partial b_q}\left[\ln\frac{b_q}{b} + \ln\frac{\sinh b}{\sinh b_q}\right] = \frac{1}{b_q} - \coth b_q$

• Then, the KL core term:

$\frac{\partial}{\partial b_q}\left[(b_q - b)S(b_q)\right] = S(b_q) + (b_q - b)\frac{dS}{db_q}$

• The linear terms:

$\frac{\partial}{\partial b_q}\left[-\sum_{j≠i}(b_j + S(b_q)J_{ij})\sigma_j\right] = -\sum_{j≠i}J_{ij}\sigma_j\frac{dS}{db_q}$

• The constants vanish.

• Now, combining all terms:

$\frac{\partial F}{\partial b_q} = \left(\frac{1}{b_q} - \coth b_q\right) + \left(S(b_q) + (b_q - b)\frac{dS}{db_q}\right) - \sum_{j≠i}J_{ij}\sigma_j\frac{dS}{db_q}$

• Substituting $S(b_q) = \coth b_q - 1/b_q$:

$= \left(\frac{1}{b_q} - \coth b_q\right) + \left(\coth b_q - \frac{1}{b_q} + (b_q - b)\frac{dS}{db_q}\right) - \sum_{j≠i}J_{ij}\sigma_j\frac{dS}{db_q}$

• Cancel terms:

$\cancel{\frac{1}{b_q}} \cancel{- \coth b_q} + \cancel{\coth b_q} \cancel{- \frac{1}{b_q}} + (b_q - b)\frac{dS}{db_q} - \sum_{j≠i}J_{ij}\sigma_j\frac{dS}{db_q}$

Gives us the final derivative:

$\frac{∂F}{∂b_q} = \left(b_q - b - \sum_{j\neq i}J_{ij}σ_j\right)\frac{dS}{db_q}$

Where $\frac{dS}{db_q} = -csch^2 b_q + \frac{1}{b_q^2}$.

Setting the derivative to zero and solving for $b_q$, we get: $b_q = b + \sum_{j \ne i} J_{ij} \sigma_j$

Now we remember that the expected value of the $\mathcal{CB}$ is the Langevin function of its bias -, $\mathbb{E}(x) = coth(b) - 1/b$ (). For simplicity, we will denote it as $L(x)$. Now we can write:
$\mathbb{E}_{q}[\sigma_i] = L(b_q) = L \left( b + \sum_{j \ne i} J_{ij} \sigma_j \right)$

Q.E.D.

Appendix 7¶

Basin geometry and adversarial robustness

We provide a geometric analysis of how attractor orthogonality and the precision parameter jointly affect adversarial robustness in the proposed framework.

Setup. Consider a network with $N$ nodes, inverse temperature $iT$, and $K$ stored attractors $\{\boldsymbol{\sigma}^{(\mu)}\}_{\mu=1}^K$ with $\|\boldsymbol{\sigma}^{(\mu)}\|^2 = N$. Under the projector weight matrix $\mathbf{J}^\dagger \approx \frac{1}{N}\sum_\mu \boldsymbol{\sigma}^{(\mu)}{\boldsymbol{\sigma}^{(\mu)}}^\top$ (zero diagonal), which the learning rule approximates, the energy at a stored attractor satisfies $E(\boldsymbol{\sigma}^{(\mu)}) \approx -N$.

Adversarial budget and representational redundancy. An adversarial perturbation $\delta\mathbf{s}$ applied to the bias vector shifts the energy difference between two attractors $\boldsymbol{\sigma}^{(\mu)}$ and $\boldsymbol{\sigma}^{(\nu)}$ by $\Delta E \approx -\delta\mathbf{s}^\top (\boldsymbol{\sigma}^{(\mu)} - \boldsymbol{\sigma}^{(\nu)})$. The most efficient attack aligns $\delta\mathbf{s}$ with the difference vector, yielding a minimum perturbation norm $\|\delta\mathbf{s}\|_{\min} \propto 1/\|\boldsymbol{\sigma}^{(\mu)} - \boldsymbol{\sigma}^{(\nu)}\|$. For orthogonal attractors ($\boldsymbol{\sigma}^{(\mu)\top}\boldsymbol{\sigma}^{(\nu)} = 0$): $\|\boldsymbol{\sigma}^{(\mu)} - \boldsymbol{\sigma}^{(\nu)}\|^2 = 2N$. For attractors with normalized overlap $m = \frac{1}{N}\boldsymbol{\sigma}^{(\mu)\top}\boldsymbol{\sigma}^{(\nu)}$: $\|\boldsymbol{\sigma}^{(\mu)} - \boldsymbol{\sigma}^{(\nu)}\|^2 = 2N(1-m)$. Since $0 \leq m < 1$, orthogonal attractors ($m = 0$) maximize the adversarial budget. Self-orthogonalization thus reduces representational redundancy (overlap between attractor states) while increasing robustness to targeted perturbations.

Structural redundancy. Although representational redundancy is minimized, structural redundancy is fully preserved: each attractor is encoded across all $O(N^2)$ synaptic weights, with each weight contributing $O(1/N)$ to the basin depth. Corrupting or removing a bounded fraction of weights therefore has a limited effect on the energy landscape, regardless of attractor orthogonality.

Prior precision as a defense. Under finite temperature, the posterior is $p(\boldsymbol{\sigma}\mid\mathbf{s}) \propto \exp\bigl\{iT\bigl[\sum_i (b_i + \delta s_i)\sigma_i + \frac{1}{2}\sum_{ij}J_{ij}\sigma_i\sigma_j\bigr]\bigr\}$. The effective basin depth scales as $iT \cdot N$, while the perturbation contributes $iT \cdot \delta\mathbf{s}^\top\boldsymbol{\sigma}$. The ratio of basin depth to perturbation magnitude grows with $iT$: high precision makes prior-encoded attractors robust to sensory-level attacks. The converse also holds: low $iT$ down-weights the prior, limiting the influence of poisoned attractors on the posterior.

Stochastic averaging. Beyond the geometric analysis, stochastic MCMC inference provides an additional defense layer. The time-averaged network response converges to the posterior mean $\mathbb{E}_{p(\boldsymbol{\sigma}\mid\mathbf{s})}[\boldsymbol{\sigma}]$, with variance decreasing as $O(1/T_{\text{eff}})$ where $T_{\text{eff}}$ is the effective number of independent samples. This averaging smooths out the effect of perturbations that shift individual samples but do not move the posterior mode across a basin boundary.

Implicit regularization against data poisoning. During learning, the anti-Hebbian term in the learning rule subtracts the network’s current prediction from the observed correlation. For a corrupted training pattern, the prediction — dominated by the clean attractor subspace — absorbs most of the pattern’s variance, leaving only a small residual update. For a poisoning fraction $\epsilon \ll 1$, the cumulative effect of poisoned updates remains bounded relative to the clean learning signal, analogous to the outlier-robustness of Bayesian posteriors under informative priors.

Appendix 8¶

Detailed comparison with classical models

We compare the proposed FEP-based attractor network (FEP-ANN) with the canonical formulations of classical Hopfield networks, Boltzmann machines, projector (pseudo-inverse) networks, and hierarchical predictive coding. Each of these model families encompasses many variants; the comparison targets the canonical baseline in each case.

Learning rule and computational complexity

• Boltzmann machines (contrastive divergence). The exact gradient of the Boltzmann log-likelihood with respect to weights $J_{ij}$ is

  $$\frac{\partial \ln p(\mathbf{v})}{\partial J_{ij}} = \langle \sigma_i \sigma_j \rangle_{\text{data}} - \langle \sigma_i \sigma_j \rangle_{\text{model}}$$

  (2)

  where $\langle \cdot \rangle_{\text{model}}$ requires sampling from the model’s equilibrium distribution — computationally intractable for large networks. Contrastive divergence (CD-$k$) approximates $\langle \cdot \rangle_{\text{model}}$ by running $k$ steps of Gibbs sampling from the data-initialized state. Each Gibbs sweep visits all $N$ nodes, and computing the input to each node requires $O(N)$ multiplications, yielding a per-step cost of $O(N^2)$ per Gibbs sweep and $O(N^2 k)$ per weight update.

• FEP-ANN. The learning rule for node $i$ uses the network’s instantaneous predicted state $\hat{\sigma}_i = L(b_i + \sum_k J_{ik} \sigma_k)$ — a single forward pass through the current weights — rather than a sample average over model-generated states. Computing $\hat{\sigma}_i$ for all nodes requires one matrix-vector product $\mathbf{J}\boldsymbol{\sigma}$ at $O(N^2)$; the outer-product updates for all weights then cost $O(N^2)$. The total per-step learning cost is therefore $O(N^2)$ — a factor of $k$ cheaper than CD-$k$. The free-running phase is eliminated entirely.

• Solenoidal flows and mixing times In a symmetric Boltzmann machine, inference proceeds via reversible Gibbs sampling; reversibility implies that probability mass can only move between attractor basins by climbing over energy barriers, leading to mixing times that scale exponentially with barrier height.
  
  In the FEP-ANN with asymmetric couplings $\mathbf{J} \neq \mathbf{J}^\top$, the dynamics decompose into a gradient (dissipative) component driven by $\mathbf{J}^\dagger = \frac{1}{2}(\mathbf{J}+\mathbf{J}^\top)$ and a solenoidal (conservative) component driven by $\mathbf{J}^- = \frac{1}{2}(\mathbf{J}-\mathbf{J}^\top)$. The solenoidal term does not alter the stationary distribution but induces persistent probability currents along iso-energy contours. These non-reversible currents transport probability mass between attractor basins without climbing energy barriers — they traverse along the iso-energy surface instead, accelerating mixing. The resulting dynamics are formally analogous to continuous normalizing flows and Markovian flow matching.

• Orthogonalization and memory capacity Classical Hopfield networks with Hebbian (outer-product) weights store $K$ patterns; cross-talk limits capacity to $K_{\max} \approx 0.14N$. Projector (pseudo-inverse) networks achieve optimal capacity $K_{\max} = N$ with error-free retrieval but require batch access and $O(N^3)$ computation. In the FEP-ANN, the Hebbian/anti-Hebbian learning rule incrementally approximates the projector matrix through online learning: the Hebbian term accumulates pattern correlations while the anti-Hebbian term subtracts the model’s current prediction, progressively decorrelating the stored representations. Memory capacity thus approaches the projector limit without batch computation or matrix inversion.

• Relationship to hierarchical predictive coding Both the FEP-ANN and hierarchical predictive coding (PC) networks minimize variational free energy via local, prediction-error-driven updates. The structural difference is topology: PC operates on a bidirectional hierarchy with directed inter-layer connections (feedforward predictions, feedback prediction errors), while the FEP-ANN uses a fully recurrent, single-layer topology with arbitrary lateral couplings. Within the deep particular partition formalism, both are instances of the same FEP-based inference architecture applied to different coupling topologies. Exploring the formal equivalences between these topologies is a promising direction for future work.

Appendix 9¶

Scaling benchmarks

We report runtime and memory capacity benchmarks for the JAX implementation across network sizes from $N = 64$ to $N = 4096$ and pattern counts from $K = 10$ to $K = N$. Training and inference wall-clock time as a function of network size $N$, confirming the expected $O(N^2)$ per-step scaling. Panel to be finalized with benchmark data from 07-simulation-scaling-jax.ipynb.

Table: Runtime benchmarks (JAX implementation, 10 updates)