Mean Field Approximation

Overview: When and Why Mean Field?

Mean field approximation is a powerful technique that appears in two seemingly different domains: statistical inference (factorizing posterior distributions) and neuroscience (replacing spike interactions with rate expectations).

One-Sentence Unification: Mean field approximation replaces many-body interactions with self-consistent expectations, whether those bodies are neurons exchanging spikes or random variables exchanging information.

The Deep Connection: Free Energy

This isn't just analogy—it's mathematical equivalence. Mean field neuroscience can be derived from maximum entropy, Gibbs measures, and statistical mechanics of spiking networks.

Rate-based mean field equations correspond to stationary points of a free energy functional
Asynchronous irregular states ≈ factorized Gibbs distributions
Active Inference literally imports variational mean field machinery into neural process theory

When Friston talks about neurons minimizing free energy, he is not using metaphor—he is exploiting the mathematical equivalence.

What You'll Learn

Topic	Key Concept	Takeaway
Single Factor POMDP	Exact Bayesian inference	No structural approximation needed when state space is small
Algorithmic vs Structural	Two kinds of approximation	Gradient descent ≠ mean field factorization
Free Energy Gradient	∂F/∂q = log q + 1 - log p	Where the +1 comes from (product rule!)
Multi Factor POMDP	Joint vs factorized inference	Correlations emerge when factors interact
Covariance Deep Dive	E[XY] - E[X]E[Y]	Why MFA forces Cov = 0
MF Derivation	Expectations over other factors	Coordinate descent from free energy
The Click Moment	Fixed points & self-consistency	Independence given the mean
Neural Mean Field	Spikes → Rates	Same idea applied to neural populations
The Isomorphism	Math equivalence	One algorithm, two languages

Single Factor POMDP: No Structural Approximation Needed

The Model

Hidden state: s ∈ { Near, Mid, Far }
Observations: Angle bucket ("Toward" 0-120°, "Sideways" 120-240°, "Away" 240-360°)
State space size: Only 3 states!

The Bayes Rule

Exact Bayes: q(s) = p(s|o) ∝ p(o|s) × p(s)

What does ∝ mean? The symbol ∝ means "proportional to". We write ∝ instead of = because the right side isn't normalized—it doesn't sum to 1. To get actual probabilities, we divide by the marginal likelihood (also called "evidence"):

p(o) = Σₛ p(o|s)×p(s)

This p(o) is the probability of the observation marginalized over all possible states. So the full equation is:

q(s) = p(o|s) × p(s) / p(o)

We often skip writing p(o) because it's just a normalizing factor that makes everything sum to 1.

Free Energy: The Variational Perspective

Instead of computing Bayes' rule directly, we can minimize free energy:

F[q] = KL(q || p(s|o)) + log p(o)

Since log p(o) is constant (it doesn't depend on q), minimizing F is equivalent to minimizing KL divergence—which means finding q that matches the true posterior!

Expanding Free Energy

F[q] = Σₛ q(s) × [ log q(s) - log p(o,s) ]
     = Σₛ q(s) × [ log q(s) - log p(o|s) - log p(s) ]

The Gradient (for descent)

∂F/∂q(s) = log q(s) + 1 - log p(o|s) - log p(s)
             ↑         ↑
          from q    product rule on q×log(q)

The +1 comes from the product rule: ∂/∂q [q×log q] = log q + 1. At equilibrium (∇F = 0), this gives us q* = p(s|o).

Try It: Bayesian Inference

Exact (Bayes Rule)Variational (Gradient Descent on F)

Near

33.3%

Mid

33.3%

Far

33.3%

Both methods give the same answer! The variational approach is "extra" here, but it's didactically essential—it shows the FEP machinery in the simplest case.

Multi Factor POMDP: When Things Get Interesting

The Model

Factor 1: Distance ∈ { Near, Mid, Far }
Factor 2: Direction ∈ { Toward, Away }
Joint state space: 3 × 2 = 6 states

Now observations depend on both factors. If you're Near and facing Toward, you get a stronger signal than if you're Near but facing Away. This creates correlationsin the posterior.

Why this matters: With only 6 joint states, we can still do exact inference. But imagine adding more factors—with 10 binary factors, you'd have 2¹⁰ = 1,024 states. With 100 factors, you'd have 2¹⁰⁰ states (more than atoms in the universe!). This is where mean field approximation becomes essential.

Exact: p(distance, direction | o) — full 6-state joint with correlations
Mean Field: q(distance) × q(direction) — only 3+2 = 5 parameters, no correlations

Try It: Compare Exact vs Mean Field

Exact (Joint)Mean Field (Factorized)

Distance Beliefs:

Near

33.3%

Mid

33.3%

Far

33.3%

Direction Beliefs:

Toward

50.0%

Away

50.0%

Full Joint (6 states):

	Toward	Away
Near	16.7%	16.7%
Mid	16.7%	16.7%
Far	16.7%	16.7%

Method: Exact (preserves correlations) Covariance: 0.0000

Notice how the marginals are similar between methods, but mean field loses the correlation information. The covariance tells us how knowing one factor informs the other—mean field discards this.

Two Kinds of Approximation

This is a crucial distinction that often causes confusion:

Type	What's Approximated	Can Recover Exact?
Algorithmic	The optimization path (gradient descent vs. closed-form)	✅ Yes (with more iterations)
Structural (Mean Field)	The posterior's form: q(s₁,s₂) = q(s₁) × q(s₂)	❌ No (correlations permanently lost)

Algorithmic Approximation (Gradient Descent)

for i in iterations:
  grad = ∂F/∂q
  q = q - learningRate × grad

This is an iterative algorithm to find the same answer as closed-form Bayes. With enough iterations, it converges to: q(s) ∝ p(o|s) × p(s)

Structural Approximation (Mean Field)

q(distance, direction) ≈ q(distance) × q(direction)

This is a restriction on the form of the solution. Even with infinite iterations, you cannot recover the true correlated posterior.

Variational Inference is the framework. Mean Field is just one choiceof tractable family: the family of fully factorized distributions q(s) = Π qᵢ(sᵢ).

The Variational Inference Hierarchy

Variational Inference (the framework)
    │
    ├── Exact VI (family = all distributions)
    │       └── Gradient descent (converges to exact)
    │       └── Closed-form Bayes rule
    │
    ├── Mean Field VI (family = factorized)
    │       └── Coordinate descent
    │
    ├── Gaussian VI (family = Gaussians)
    │
    └── ... other restricted families

In our single-factor POMDP, gradient descent is "extra" computationally—but it'sessential didactically because it shows the Free Energy Principle in action: an agent literally minimizing F!

Covariance Deep Dive: Why Mean Field Forces Cov = 0

Understanding covariance is essential for grasping what mean field loses.

Intuitive Definitions

Concept	Intuitive Name
Variance	"Average squared wiggle of X from its mean"
Standard Deviation	"Average wiggle of X from its mean" (standardized)
Covariance	"Average co-wiggle of X and Y from their means"
Correlation	"Standardized co-wiggle" (between -1 and +1)

The Covariance Formula

Cov(X, Y) = E[XY] - E[X]×E[Y]

This asks: "How much does E[XY] differ from what it WOULD be if X and Y were independent?"

When X is...	And Y is...	Product is...	Meaning
Above mean (+)	Above mean (+)	Positive	They move together!
Above mean (+)	Below mean (-)	Negative	They move opposite!
Below mean (-)	Below mean (-)	Positive	They move together!
Below mean (-)	Above mean (+)	Negative	They move opposite!

Cov > 0: Products mostly positive → they move together
Cov < 0: Products mostly negative → they move opposite
Cov ≈ 0: Positives and negatives cancel → independent

Why Mean Field Forces Covariance = 0

Here's the key insight. Mean field says: q(d, dir) = q(d) × q(dir)

This defines the joint as a product of marginals. By construction:

E[distance × direction] 
= Σ_d Σ_dir  d × dir × q(d, dir)
= Σ_d Σ_dir  d × dir × q(d) × q(dir)  [mean field!]
= (Σ_d d × q(d)) × (Σ_dir dir × q(dir))
= E[d] × E[dir]

Therefore:

Cov = E[d × dir] - E[d] × E[dir]
    = E[d] × E[dir] - E[d] × E[dir]
    = 0  ← ALWAYS zero in mean field!

Mean field doesn't "discover" that variables are independent—it assumes they are. The factorization forces E[XY] = E[X]×E[Y], which is the definition of independence!

The Correlation Bound (Cauchy-Schwarz)

Here's something beautiful: |Cov(X,Y)| ≤ StdDev(X) × StdDev(Y). Always!

Dividing by this bound gives correlation: a unitless value between -1 and +1.

Correlation = Cov(X,Y) / (σₓ × σᵧ)

This is the "standardized" covariance—comparable across different variable pairs.

Mean Field Derivation: From Free Energy to Coordinate Descent

Let's derive the mean field update equations from the variational free energy. This shows exactly where the "expectation over other factors" comes from.

Setup: Two Hidden Factors

d   ∈ {0, 1, 2}     (Near, Mid, Far)
dir ∈ {0, 1}        (Toward, Away)
Joint state space: 6 combinations (d, dir)
Observation: o (which angle bucket we see)
Generative model: p(o, d, dir) = p(o | d, dir) × p(d, dir)

Step 1: Write the Free Energy

F[q] = Σ_d Σ_dir q(d, dir) × [ log q(d, dir) - log p(o, d, dir) ]

Step 2: Apply Mean Field Assumption

The assumption: q(d, dir) = q_d(d) × q_dir(dir)

Substitute and use log product rule:

F = Σ_d Σ_dir [q_d(d) × q_dir(dir)] × [ log q_d(d) + log q_dir(dir) - log p(o, d, dir) ]

Step 3: Expand and Simplify

F = Σ_d q_d(d) log q_d(d)      [sum over dir = 1]
  + Σ_dir q_dir(dir) log q_dir(dir)   [sum over d = 1]
  - Σ_d Σ_dir q_d(d) q_dir(dir) log p(o, d, dir)

Step 4: Take Derivative w.r.t. q_d(d)

∂F/∂q_d(d) = log q_d(d) + 1 - Σ_dir q_dir(dir) log p(o, d, dir)
                              └────────────────────────────────┘
                                   This is an EXPECTATION!

Step 5: Set Gradient to Zero

At the minimum:

log q_d(d) = E_q(dir) [log p(o, d, dir)] - 1 + const

The -1 + const is absorbed in normalization, so:

q_d(d) ∝ exp( E_q(dir) [log p(o, d, dir)] )

The Mean Field Algorithm

1. Initialize:
   q_d = [1/3, 1/3, 1/3]
   q_dir = [1/2, 1/2]

2. Repeat until converged:

   a. Update q_d (holding q_dir fixed):
      For each d:
        log q_d(d) = Σ_dir q_dir(dir) × log p(o, d, dir)
      Normalize: q_d = softmax(log q_d)

   b. Update q_dir (holding q_d fixed):
      For each dir:
        log q_dir(dir) = Σ_d q_d(d) × log p(o, d, dir)
      Normalize: q_dir = softmax(log q_dir)

3. Return q_d, q_dir

Compare: Single vs Multi-Factor Gradient

Single Factor	Multi-Factor Mean Field
∂F/∂q(s) = log q(s) + 1 - log p(o\|s) - log p(s)	∂F/∂q_d(d) = log q_d(d) + 1 - E_q(dir)[log p(o,d,dir)]

The only difference: Instead of the actual log p(...), we use theexpected log p(...) over our current belief about the other factor! This is the mean field move: replace interaction with expectation.

Mean Field Approximation: Factorizing the Posterior

Mean Field Assumption: q(distance, direction) ≈ q(distance) × q(direction)

Mean field approximation replaces the joint posterior with a product of marginals. This makes inference tractable but forces correlations to zero.

The Algorithm: Coordinate Descent

Initialize q(distance) and q(direction) uniformly
Update q(distance) using expectation over q(direction)
Update q(direction) using expectation over q(distance)
Repeat until self-consistent (fixed point)

Compare: Exact vs Mean Field (Side by Side)

📊 Exact Marginals

Near

33%

Mid

33%

Far

33%

Toward

50%

Away

50%

Cov: 0.0000

📊 Mean Field

Near

33%

Mid

33%

Far

33%

Toward

50%

Away

50%

Cov: 0.0000 (always!)

Mean field gets the marginals approximately right but completely discards correlations. The factorization q(d)×q(dir) forces Cov(d, dir) = 0 by construction.

La Camera's Binary Neurons: A Neural Model

Before we see "the click" where neural and statistical mean field become one, let's introduce the neural model from La Camera (2022).

Binary Neuron Network

Each neuron i has a binary state xᵢ ∈ { 0, 1 } (silent or spiking). Neurons are connected by synaptic weights wᵢⱼ.

Variable Reference

Symbol	Meaning	Domain
xᵢ(t)	Spike state of neuron i at time t	{ 0, 1 }
rᵢ	Firing rate of neuron i (time-averaged)	[0, 1]
wᵢⱼ	Synaptic weight from neuron j to i	ℝ
Iᵢ	Total input current to neuron i	ℝ
I_ext	External input (sensory drive)	ℝ
σ(·)	Sigmoid transfer function	(0, 1)
β	Inverse temperature (sharpness)	ℝ⁺
N	Number of neurons in the network	ℕ

The Exact Dynamics

Input:      Iᵢ(t) = Σⱼ wᵢⱼ × xⱼ(t) + I_ext
Spike prob: P(xᵢ = 1) = σ(Iᵢ) = 1 / (1 + exp(-2β × Iᵢ))
Update:     xᵢ(t+1) = 1 with probability P(xᵢ = 1)

Each neuron's firing depends on who actually spiked at time t. This creates complex dependencies and correlations.

The Mean Field Approximation

Replace spikes with rates:
  Iᵢ ≈ Σⱼ wᵢⱼ × rⱼ + I_ext

Self-consistency equation:
  rᵢ = σ(Σⱼ wᵢⱼ × rⱼ + I_ext)

Or in vector form:
  r = σ(W × r + I_ext)

Instead of tracking binary spikes xⱼ(t), we use expected rates rⱼ = ⟨xⱼ⟩. The input becomes deterministic—no more fluctuations from individual spikes!

What we lose: Spike correlations! If neurons tend to fire together (synchrony), mean field misses this. It only sees the average activity.

Why This Matters

Exact simulation: Track 2^N possible spike patterns (intractable for N > 30)
Mean field: Track N firing rates (trivially tractable for any N)
Trade-off: Speed for accuracy—correlations are discarded

This is the same move as variational mean field! Replace interactions (spikes xⱼ) with expectations (rates rⱼ), then solve for self-consistent values.

The Click: Self-Consistency and Fixed Points

This is where mean field stops being a formal trick and starts to feel like a controlled lie. Let's understand the deep structure.

The Core Problem (All Three Domains)

You have many interacting things. Tracking all interactions exactly is impossible.

Domain	The "Things"	The "Interactions"	Why Intractable
Statistics	Random variables s₁, s₂, ...	Joint dependencies	Exponential state space
Neuroscience	Neurons 1, 2, ..., N	Synaptic connections wᵢⱼ	2^N spike patterns
Physics	Particles, spins	Forces, couplings	10²³ atoms

The Mean Field Move (Same in All!)

Replace interactions with averages.

Domain	Exact	Mean Field
Statistics	q(s₁,s₂) - full joint	q(s₁) × q(s₂) - product
Neuroscience	I = Σ wᵢⱼ × xⱼ(t) - spikes	I = Σ wᵢⱼ × rⱼ - rates
Physics	Each atom feels all others	Each feels average field

What is "Self-Consistency"?

Think of it as a feedback loop:

┌─────────────────────────────────────────┐
│                                         │
│   YOUR BELIEF  ──────►  WHAT YOU SEE    │
│       ↑                      │          │
│       │                      │          │
│       └──────────────────────┘          │
│         (must agree!)                   │
│                                         │
└─────────────────────────────────────────┘

In Neurons:

Firing rate r_i ──► determines input I_j = Σ w_ji × r_i
                          │
                          ▼
                     determines firing rate r_j = σ(I_j)
                          │
                          ▼
                     determines input I_i = Σ w_ij × r_j  
                          │
                          ▼
                     back to r_i!

Self-consistent = the loop closes. What goes in equals what comes out. This is a fixed point.

In Variational Inference:

Your belief q(s₁) ──► determines E[s₁]
                          │
                          ▼
                     affects how you update q(s₂)
                          │
                          ▼
                     which affects E[s₂]
                          │
                          ▼
                     which affects how you update q(s₁)
                          │
                          ▼
                     back to the start!

The Key Insight: Independence Given the Mean

This is perfectly stated:

Independence given the mean: Mean field assumes that given the population rate (or expectation), individual neurons (or variables) fluctuate independently. We represent only the marginal tendencies and let expectations stand in for interactions.

Dynamical Systems Language

Term	Meaning	Example
Fixed point	A state that doesn't change	r* such that r* = σ(Wr* + I)
Iteration	Repeat until stable	r(t+1) = σ(Wr(t) + I)
Convergence	Getting closer to fixed point	r(t) → r* as t → ∞
Self-consistent	Fixed point exists & meaningful	The physics/stats makes sense

The One-Sentence Summary

Mean field theory replaces microscopic interactions with macroscopic averages, then demands that those averages be self-consistent with the behavior they produce.

This sentence works for variational Bayes, spiking neural networks, magnetic materials, economics (mean field games), epidemiology (SIR models)—it's the same mathematical structure everywhere.

Neural Mean Field: Spikes → Rates

The Key Move

Exact: I_i = Σ_j w_ij × x_j(t) — depends on actual spikes
Mean Field: I_i ≈ Σ_j w_ij × ⟨r_j⟩ — depends on expected rates

In a network of N neurons, each neuron i receives input from all others. The exact input depends on who actually spiked. Mean field replaces this with the expected input based on firing rates.

The Self-Consistency Equation

Exact simulation equations:
  I_i(t) = Σ_{j≠i} w_ij × x_j(t) + I_ext
  P(spike) = σ(I_i) = 1 / (1 + exp(-2β × I_i))
  x_i(t+1) = 1 with probability P(spike)

Mean field self-consistency:
  r = σ(W × r + I_ext)

For a HOMOGENEOUS network (all neurons identical):
  r* = σ(J × r* + I_ext)

where J = N × w̄ is the TOTAL EFFECTIVE COUPLING
      N = number of neurons
      w̄ = average synaptic weight per connection

Each neuron's rate depends on others' rates. Iterate until consistent.

Step-by-Step Iteration Trace

Watch mean field converge to its fixed point

Compare: Exact Simulation vs Mean Field

100 neurons, 500 time steps

Click "Run Comparison" to see exact vs mean field...

A Deep Insight: Attractors Are Nowhere in the Equations!

The equation r = σ(J × r + I_ext) has r on both sides. It doesn't guarantee a fixed point exists—you have to prove it.

Brouwer's Fixed Point Theorem: Any continuous function from a closed bounded set to itself has at least one fixed point. Since σ: ℝ → (0,1) is continuous and maps [0,1] → (0,1) ⊂ [0,1], a fixed point must exist.

Bistability: When There Are Multiple Fixed Points

For strong enough coupling J, the sigmoid curve can intersect the diagonal three times—givingtwo stable attractors and one unstable fixed point in between!

          r
          ↑
     1.0  │           .---------  ← σ(Jr + I)
          │         .´
          │       .´
          │     .´    *  ← Fixed point (stable)
          │    /
          │   /
          │  /   *  ← Fixed point (unstable)  
          │ /
          │/
          │    *  ← Fixed point (stable)
     0.0  └────────────────────→ r
          0                    1

Which attractor do you find? It depends on initial conditions! This isbistability—the system has memory. In the brain, bistability can represent working memory, decision making, or attention.

Mean field predicts qualitative phase transitions: weak coupling → one fixed point, strong coupling → bistability. These predictions are often correct even when quantitative details are wrong!

The Isomorphism: One Algorithm, Two Languages

The mathematical structure of mean field is identical in both domains:

Concept	Neuroscience	Variational Inference
Individual unit in the system	Neuron i	Variable sᵢ
Discrete state of each unit	Spike xᵢ ∈ {0,1}	State sᵢ ∈ {1,...,K}
Expected value / distribution	Firing rate rᵢ = ⟨xᵢ⟩	Marginal qᵢ(sᵢ)
Interaction strength	Weight wᵢⱼ	Coupling log p(sᵢ\|sⱼ)
Aggregated influence from others	Input Iᵢ = Σⱼ wᵢⱼrⱼ	log qᵢ = E_{qⱼ}[log p]
Nonlinear activation	Transfer φ(I) = sigmoid	Transfer = softmax
What mean field ignores	Noise η	Discarded correlations
Independence assumption	Asynchronous state	Factorized q = Πᵢ qᵢ
Self-consistency equation	r = φ(Wr + Iₑₓₜ)	Fixed point of coord descent

The Algorithm (Same in Both!)

Step	Neural	Variational
Initialize Expectations	r_i ← 0.5 (random rates)	q_i ← uniform (uninformative marginals)
Compute Input	I_i = Σ_j w_ij × r_j + I_ext	log q_i = E_{q_{-i}}[log p(s_i, s_{-i})]
Apply Transfer	r_i ← sigmoid(I_i)	q_i ← softmax(log q_i)
Check Convergence	max\|r_new - r_old\| < ε	max\|q_new - q_old\| < ε

Key Insights

The Core Idea: Replace INTERACTIONS with EXPECTATIONS, then demand SELF-CONSISTENCY.

What We Lose: Correlations, synchrony, and joint patterns are discarded. Covariance is forced to zero.

What We Gain: Tractability! Instead of exponential complexity, we get polynomial or linear updates.

When It Works: Best when correlations are weak, populations are large, or we only care about means.

The Free Energy Connection: Both are minimizing a free energy functional. Variational F[q] = KL(q||p), Neural F = energy - entropy.

Quick Reference

When to Use Mean Field

Scenario	Mean Field Needed?	Why?
Single factor, small state space	No	Exact inference is tractable
Multiple interacting factors	Often	Joint space grows exponentially
Neural population with correlations	Yes	Spike-by-spike tracking intractable
Weak correlations, care about means	Good fit	Mean field error is small
Strong correlations matter	Problematic	Mean field discards what you need

Key Equations

Concept	Equation
Bayes Rule	`p(s\|o) ∝ p(o\|s) × p(s)`
Free Energy	`F[q] = Σ q(s)[log q(s) - log p(o,s)]`
Free Energy Gradient	`∂F/∂q(s) = log q(s) + 1 - log p(o\|s) - log p(s)`
Mean Field Factorization	`q(s₁,s₂) ≈ q(s₁) × q(s₂)`
Neural Self-Consistency	`r = σ(W × r + I_ext)`
Covariance	`Cov(X,Y) = E[XY] - E[X]×E[Y]`

The One-Liner

Mean field approximation replaces many-body interactions with self-consistent expectations, whether those bodies are neurons exchanging spikes or random variables exchanging information.

Loading tutorial…

Overview: When and Why Mean Field?

The Deep Connection: Free Energy

This isn't just analogy—it's mathematical equivalence. Mean field neuroscience can be derived from maximum entropy, Gibbs measures, and statistical mechanics of spiking networks.

Rate-based mean field equations correspond to stationary points of a free energy functional
Asynchronous irregular states ≈ factorized Gibbs distributions
Active Inference literally imports variational mean field machinery into neural process theory

When Friston talks about neurons minimizing free energy, he is not using metaphor—he is exploiting the mathematical equivalence.

What You'll Learn

Topic	Key Concept	Takeaway
Single Factor POMDP	Exact Bayesian inference	No structural approximation needed when state space is small
Algorithmic vs Structural	Two kinds of approximation	Gradient descent ≠ mean field factorization
Free Energy Gradient	∂F/∂q = log q + 1 - log p	Where the +1 comes from (product rule!)
Multi Factor POMDP	Joint vs factorized inference	Correlations emerge when factors interact
Covariance Deep Dive	E[XY] - E[X]E[Y]	Why MFA forces Cov = 0
MF Derivation	Expectations over other factors	Coordinate descent from free energy
The Click Moment	Fixed points & self-consistency	Independence given the mean
Neural Mean Field	Spikes → Rates	Same idea applied to neural populations
The Isomorphism	Math equivalence	One algorithm, two languages

Single Factor POMDP: No Structural Approximation Needed

The Model

Hidden state: s ∈ { Near, Mid, Far }
Observations: Angle bucket ("Toward" 0-120°, "Sideways" 120-240°, "Away" 240-360°)
State space size: Only 3 states!

The Bayes Rule

Exact Bayes: q(s) = p(s|o) ∝ p(o|s) × p(s)

Free Energy: The Variational Perspective

Instead of computing Bayes' rule directly, we can minimize free energy:

F[q] = KL(q || p(s|o)) + log p(o)

Since log p(o) is constant (it doesn't depend on q), minimizing F is equivalent to minimizing KL divergence—which means finding q that matches the true posterior!

Expanding Free Energy

F[q] = Σₛ q(s) × [ log q(s) - log p(o,s) ]
     = Σₛ q(s) × [ log q(s) - log p(o|s) - log p(s) ]

The Gradient (for descent)

∂F/∂q(s) = log q(s) + 1 - log p(o|s) - log p(s)
             ↑         ↑
          from q    product rule on q×log(q)

The +1 comes from the product rule: ∂/∂q [q×log q] = log q + 1. At equilibrium (∇F = 0), this gives us q* = p(s|o).

Try It: Bayesian Inference

Exact (Bayes Rule)Variational (Gradient Descent on F)

Near

33.3%

Mid

33.3%

Far

33.3%

Both methods give the same answer! The variational approach is "extra" here, but it's didactically essential—it shows the FEP machinery in the simplest case.

Multi Factor POMDP: When Things Get Interesting

The Model

Factor 1: Distance ∈ { Near, Mid, Far }
Factor 2: Direction ∈ { Toward, Away }
Joint state space: 3 × 2 = 6 states

Now observations depend on both factors. If you're Near and facing Toward, you get a stronger signal than if you're Near but facing Away. This creates correlationsin the posterior.

Exact: p(distance, direction | o) — full 6-state joint with correlations
Mean Field: q(distance) × q(direction) — only 3+2 = 5 parameters, no correlations

Try It: Compare Exact vs Mean Field

Exact (Joint)Mean Field (Factorized)

Distance Beliefs:

Near

33.3%

Mid

33.3%

Far

33.3%

Direction Beliefs:

Toward

50.0%

Away

50.0%

Full Joint (6 states):

	Toward	Away
Near	16.7%	16.7%
Mid	16.7%	16.7%
Far	16.7%	16.7%

Method: Exact (preserves correlations) Covariance: 0.0000

Two Kinds of Approximation

This is a crucial distinction that often causes confusion:

Type	What's Approximated	Can Recover Exact?
Algorithmic	The optimization path (gradient descent vs. closed-form)	✅ Yes (with more iterations)
Structural (Mean Field)	The posterior's form: q(s₁,s₂) = q(s₁) × q(s₂)	❌ No (correlations permanently lost)

Algorithmic Approximation (Gradient Descent)

for i in iterations:
  grad = ∂F/∂q
  q = q - learningRate × grad

This is an iterative algorithm to find the same answer as closed-form Bayes. With enough iterations, it converges to: q(s) ∝ p(o|s) × p(s)

Structural Approximation (Mean Field)

q(distance, direction) ≈ q(distance) × q(direction)

This is a restriction on the form of the solution. Even with infinite iterations, you cannot recover the true correlated posterior.

Variational Inference is the framework. Mean Field is just one choiceof tractable family: the family of fully factorized distributions q(s) = Π qᵢ(sᵢ).

The Variational Inference Hierarchy

Variational Inference (the framework)
    │
    ├── Exact VI (family = all distributions)
    │       └── Gradient descent (converges to exact)
    │       └── Closed-form Bayes rule
    │
    ├── Mean Field VI (family = factorized)
    │       └── Coordinate descent
    │
    ├── Gaussian VI (family = Gaussians)
    │
    └── ... other restricted families

In our single-factor POMDP, gradient descent is "extra" computationally—but it'sessential didactically because it shows the Free Energy Principle in action: an agent literally minimizing F!

Covariance Deep Dive: Why Mean Field Forces Cov = 0

Understanding covariance is essential for grasping what mean field loses.

Intuitive Definitions

Concept	Intuitive Name
Variance	"Average squared wiggle of X from its mean"
Standard Deviation	"Average wiggle of X from its mean" (standardized)
Covariance	"Average co-wiggle of X and Y from their means"
Correlation	"Standardized co-wiggle" (between -1 and +1)

The Covariance Formula

Cov(X, Y) = E[XY] - E[X]×E[Y]

This asks: "How much does E[XY] differ from what it WOULD be if X and Y were independent?"

When X is...	And Y is...	Product is...	Meaning
Above mean (+)	Above mean (+)	Positive	They move together!
Above mean (+)	Below mean (-)	Negative	They move opposite!
Below mean (-)	Below mean (-)	Positive	They move together!
Below mean (-)	Above mean (+)	Negative	They move opposite!

Cov > 0: Products mostly positive → they move together
Cov < 0: Products mostly negative → they move opposite
Cov ≈ 0: Positives and negatives cancel → independent

Why Mean Field Forces Covariance = 0

Here's the key insight. Mean field says: q(d, dir) = q(d) × q(dir)

This defines the joint as a product of marginals. By construction:

E[distance × direction] 
= Σ_d Σ_dir  d × dir × q(d, dir)
= Σ_d Σ_dir  d × dir × q(d) × q(dir)  [mean field!]
= (Σ_d d × q(d)) × (Σ_dir dir × q(dir))
= E[d] × E[dir]

Therefore:

Cov = E[d × dir] - E[d] × E[dir]
    = E[d] × E[dir] - E[d] × E[dir]
    = 0  ← ALWAYS zero in mean field!

Mean field doesn't "discover" that variables are independent—it assumes they are. The factorization forces E[XY] = E[X]×E[Y], which is the definition of independence!

The Correlation Bound (Cauchy-Schwarz)

Here's something beautiful: |Cov(X,Y)| ≤ StdDev(X) × StdDev(Y). Always!

Dividing by this bound gives correlation: a unitless value between -1 and +1.

Correlation = Cov(X,Y) / (σₓ × σᵧ)

This is the "standardized" covariance—comparable across different variable pairs.

Mean Field Derivation: From Free Energy to Coordinate Descent

Let's derive the mean field update equations from the variational free energy. This shows exactly where the "expectation over other factors" comes from.

Setup: Two Hidden Factors

d   ∈ {0, 1, 2}     (Near, Mid, Far)
dir ∈ {0, 1}        (Toward, Away)
Joint state space: 6 combinations (d, dir)
Observation: o (which angle bucket we see)
Generative model: p(o, d, dir) = p(o | d, dir) × p(d, dir)

Step 1: Write the Free Energy

F[q] = Σ_d Σ_dir q(d, dir) × [ log q(d, dir) - log p(o, d, dir) ]

Step 2: Apply Mean Field Assumption

The assumption: q(d, dir) = q_d(d) × q_dir(dir)

Substitute and use log product rule:

F = Σ_d Σ_dir [q_d(d) × q_dir(dir)] × [ log q_d(d) + log q_dir(dir) - log p(o, d, dir) ]

Step 3: Expand and Simplify

F = Σ_d q_d(d) log q_d(d)      [sum over dir = 1]
  + Σ_dir q_dir(dir) log q_dir(dir)   [sum over d = 1]
  - Σ_d Σ_dir q_d(d) q_dir(dir) log p(o, d, dir)

Step 4: Take Derivative w.r.t. q_d(d)

∂F/∂q_d(d) = log q_d(d) + 1 - Σ_dir q_dir(dir) log p(o, d, dir)
                              └────────────────────────────────┘
                                   This is an EXPECTATION!

Step 5: Set Gradient to Zero

At the minimum:

log q_d(d) = E_q(dir) [log p(o, d, dir)] - 1 + const

The -1 + const is absorbed in normalization, so:

q_d(d) ∝ exp( E_q(dir) [log p(o, d, dir)] )

The Mean Field Algorithm

1. Initialize:
   q_d = [1/3, 1/3, 1/3]
   q_dir = [1/2, 1/2]

2. Repeat until converged:

   a. Update q_d (holding q_dir fixed):
      For each d:
        log q_d(d) = Σ_dir q_dir(dir) × log p(o, d, dir)
      Normalize: q_d = softmax(log q_d)

   b. Update q_dir (holding q_d fixed):
      For each dir:
        log q_dir(dir) = Σ_d q_d(d) × log p(o, d, dir)
      Normalize: q_dir = softmax(log q_dir)

3. Return q_d, q_dir

Compare: Single vs Multi-Factor Gradient

Single Factor	Multi-Factor Mean Field
∂F/∂q(s) = log q(s) + 1 - log p(o\|s) - log p(s)	∂F/∂q_d(d) = log q_d(d) + 1 - E_q(dir)[log p(o,d,dir)]

Mean Field Approximation: Factorizing the Posterior

Mean Field Assumption: q(distance, direction) ≈ q(distance) × q(direction)

Mean field approximation replaces the joint posterior with a product of marginals. This makes inference tractable but forces correlations to zero.

The Algorithm: Coordinate Descent

Initialize q(distance) and q(direction) uniformly
Update q(distance) using expectation over q(direction)
Update q(direction) using expectation over q(distance)
Repeat until self-consistent (fixed point)

Compare: Exact vs Mean Field (Side by Side)

📊 Exact Marginals

Near

33%

Mid

33%

Far

33%

Toward

50%

Away

50%

Cov: 0.0000

📊 Mean Field

Near

33%

Mid

33%

Far

33%

Toward

50%

Away

50%

Cov: 0.0000 (always!)

Mean field gets the marginals approximately right but completely discards correlations. The factorization q(d)×q(dir) forces Cov(d, dir) = 0 by construction.

La Camera's Binary Neurons: A Neural Model

Before we see "the click" where neural and statistical mean field become one, let's introduce the neural model from La Camera (2022).

Binary Neuron Network

Each neuron i has a binary state xᵢ ∈ { 0, 1 } (silent or spiking). Neurons are connected by synaptic weights wᵢⱼ.

Variable Reference

Symbol	Meaning	Domain
xᵢ(t)	Spike state of neuron i at time t	{ 0, 1 }
rᵢ	Firing rate of neuron i (time-averaged)	[0, 1]
wᵢⱼ	Synaptic weight from neuron j to i	ℝ
Iᵢ	Total input current to neuron i	ℝ
I_ext	External input (sensory drive)	ℝ
σ(·)	Sigmoid transfer function	(0, 1)
β	Inverse temperature (sharpness)	ℝ⁺
N	Number of neurons in the network	ℕ

The Exact Dynamics

Input:      Iᵢ(t) = Σⱼ wᵢⱼ × xⱼ(t) + I_ext
Spike prob: P(xᵢ = 1) = σ(Iᵢ) = 1 / (1 + exp(-2β × Iᵢ))
Update:     xᵢ(t+1) = 1 with probability P(xᵢ = 1)

Each neuron's firing depends on who actually spiked at time t. This creates complex dependencies and correlations.

The Mean Field Approximation

Replace spikes with rates:
  Iᵢ ≈ Σⱼ wᵢⱼ × rⱼ + I_ext

Self-consistency equation:
  rᵢ = σ(Σⱼ wᵢⱼ × rⱼ + I_ext)

Or in vector form:
  r = σ(W × r + I_ext)

Instead of tracking binary spikes xⱼ(t), we use expected rates rⱼ = ⟨xⱼ⟩. The input becomes deterministic—no more fluctuations from individual spikes!

What we lose: Spike correlations! If neurons tend to fire together (synchrony), mean field misses this. It only sees the average activity.

Why This Matters

Exact simulation: Track 2^N possible spike patterns (intractable for N > 30)
Mean field: Track N firing rates (trivially tractable for any N)
Trade-off: Speed for accuracy—correlations are discarded

This is the same move as variational mean field! Replace interactions (spikes xⱼ) with expectations (rates rⱼ), then solve for self-consistent values.

The Click: Self-Consistency and Fixed Points

This is where mean field stops being a formal trick and starts to feel like a controlled lie. Let's understand the deep structure.

The Core Problem (All Three Domains)

You have many interacting things. Tracking all interactions exactly is impossible.

Domain	The "Things"	The "Interactions"	Why Intractable
Statistics	Random variables s₁, s₂, ...	Joint dependencies	Exponential state space
Neuroscience	Neurons 1, 2, ..., N	Synaptic connections wᵢⱼ	2^N spike patterns
Physics	Particles, spins	Forces, couplings	10²³ atoms

The Mean Field Move (Same in All!)

Replace interactions with averages.

Domain	Exact	Mean Field
Statistics	q(s₁,s₂) - full joint	q(s₁) × q(s₂) - product
Neuroscience	I = Σ wᵢⱼ × xⱼ(t) - spikes	I = Σ wᵢⱼ × rⱼ - rates
Physics	Each atom feels all others	Each feels average field

What is "Self-Consistency"?

Think of it as a feedback loop:

┌─────────────────────────────────────────┐
│                                         │
│   YOUR BELIEF  ──────►  WHAT YOU SEE    │
│       ↑                      │          │
│       │                      │          │
│       └──────────────────────┘          │
│         (must agree!)                   │
│                                         │
└─────────────────────────────────────────┘

In Neurons:

Firing rate r_i ──► determines input I_j = Σ w_ji × r_i
                          │
                          ▼
                     determines firing rate r_j = σ(I_j)
                          │
                          ▼
                     determines input I_i = Σ w_ij × r_j  
                          │
                          ▼
                     back to r_i!

Self-consistent = the loop closes. What goes in equals what comes out. This is a fixed point.

In Variational Inference:

Your belief q(s₁) ──► determines E[s₁]
                          │
                          ▼
                     affects how you update q(s₂)
                          │
                          ▼
                     which affects E[s₂]
                          │
                          ▼
                     which affects how you update q(s₁)
                          │
                          ▼
                     back to the start!

The Key Insight: Independence Given the Mean

This is perfectly stated:

Dynamical Systems Language

Term	Meaning	Example
Fixed point	A state that doesn't change	r* such that r* = σ(Wr* + I)
Iteration	Repeat until stable	r(t+1) = σ(Wr(t) + I)
Convergence	Getting closer to fixed point	r(t) → r* as t → ∞
Self-consistent	Fixed point exists & meaningful	The physics/stats makes sense

The One-Sentence Summary

Mean field theory replaces microscopic interactions with macroscopic averages, then demands that those averages be self-consistent with the behavior they produce.

This sentence works for variational Bayes, spiking neural networks, magnetic materials, economics (mean field games), epidemiology (SIR models)—it's the same mathematical structure everywhere.

Neural Mean Field: Spikes → Rates

The Key Move

Exact: I_i = Σ_j w_ij × x_j(t) — depends on actual spikes
Mean Field: I_i ≈ Σ_j w_ij × ⟨r_j⟩ — depends on expected rates

In a network of N neurons, each neuron i receives input from all others. The exact input depends on who actually spiked. Mean field replaces this with the expected input based on firing rates.

The Self-Consistency Equation

Exact simulation equations:
  I_i(t) = Σ_{j≠i} w_ij × x_j(t) + I_ext
  P(spike) = σ(I_i) = 1 / (1 + exp(-2β × I_i))
  x_i(t+1) = 1 with probability P(spike)

Mean field self-consistency:
  r = σ(W × r + I_ext)

For a HOMOGENEOUS network (all neurons identical):
  r* = σ(J × r* + I_ext)

where J = N × w̄ is the TOTAL EFFECTIVE COUPLING
      N = number of neurons
      w̄ = average synaptic weight per connection

Each neuron's rate depends on others' rates. Iterate until consistent.

Step-by-Step Iteration Trace

Watch mean field converge to its fixed point

Compare: Exact Simulation vs Mean Field

100 neurons, 500 time steps

Click "Run Comparison" to see exact vs mean field...

A Deep Insight: Attractors Are Nowhere in the Equations!

The equation r = σ(J × r + I_ext) has r on both sides. It doesn't guarantee a fixed point exists—you have to prove it.

Bistability: When There Are Multiple Fixed Points

For strong enough coupling J, the sigmoid curve can intersect the diagonal three times—givingtwo stable attractors and one unstable fixed point in between!

          r
          ↑
     1.0  │           .---------  ← σ(Jr + I)
          │         .´
          │       .´
          │     .´    *  ← Fixed point (stable)
          │    /
          │   /
          │  /   *  ← Fixed point (unstable)  
          │ /
          │/
          │    *  ← Fixed point (stable)
     0.0  └────────────────────→ r
          0                    1

Which attractor do you find? It depends on initial conditions! This isbistability—the system has memory. In the brain, bistability can represent working memory, decision making, or attention.

The Isomorphism: One Algorithm, Two Languages

The mathematical structure of mean field is identical in both domains:

Concept	Neuroscience	Variational Inference
Individual unit in the system	Neuron i	Variable sᵢ
Discrete state of each unit	Spike xᵢ ∈ {0,1}	State sᵢ ∈ {1,...,K}
Expected value / distribution	Firing rate rᵢ = ⟨xᵢ⟩	Marginal qᵢ(sᵢ)
Interaction strength	Weight wᵢⱼ	Coupling log p(sᵢ\|sⱼ)
Aggregated influence from others	Input Iᵢ = Σⱼ wᵢⱼrⱼ	log qᵢ = E_{qⱼ}[log p]
Nonlinear activation	Transfer φ(I) = sigmoid	Transfer = softmax
What mean field ignores	Noise η	Discarded correlations
Independence assumption	Asynchronous state	Factorized q = Πᵢ qᵢ
Self-consistency equation	r = φ(Wr + Iₑₓₜ)	Fixed point of coord descent

The Algorithm (Same in Both!)

Step	Neural	Variational
Initialize Expectations	r_i ← 0.5 (random rates)	q_i ← uniform (uninformative marginals)
Compute Input	I_i = Σ_j w_ij × r_j + I_ext	log q_i = E_{q_{-i}}[log p(s_i, s_{-i})]
Apply Transfer	r_i ← sigmoid(I_i)	q_i ← softmax(log q_i)
Check Convergence	max\|r_new - r_old\| < ε	max\|q_new - q_old\| < ε

Key Insights

The Core Idea: Replace INTERACTIONS with EXPECTATIONS, then demand SELF-CONSISTENCY.

What We Lose: Correlations, synchrony, and joint patterns are discarded. Covariance is forced to zero.

What We Gain: Tractability! Instead of exponential complexity, we get polynomial or linear updates.

When It Works: Best when correlations are weak, populations are large, or we only care about means.

The Free Energy Connection: Both are minimizing a free energy functional. Variational F[q] = KL(q||p), Neural F = energy - entropy.

Quick Reference

When to Use Mean Field

Scenario	Mean Field Needed?	Why?
Single factor, small state space	No	Exact inference is tractable
Multiple interacting factors	Often	Joint space grows exponentially
Neural population with correlations	Yes	Spike-by-spike tracking intractable
Weak correlations, care about means	Good fit	Mean field error is small
Strong correlations matter	Problematic	Mean field discards what you need

Key Equations

Concept	Equation
Bayes Rule	`p(s\|o) ∝ p(o\|s) × p(s)`
Free Energy	`F[q] = Σ q(s)[log q(s) - log p(o,s)]`
Free Energy Gradient	`∂F/∂q(s) = log q(s) + 1 - log p(o\|s) - log p(s)`
Mean Field Factorization	`q(s₁,s₂) ≈ q(s₁) × q(s₂)`
Neural Self-Consistency	`r = σ(W × r + I_ext)`
Covariance	`Cov(X,Y) = E[XY] - E[X]×E[Y]`

The One-Liner

Mean field approximation replaces many-body interactions with self-consistent expectations, whether those bodies are neurons exchanging spikes or random variables exchanging information.