### New paper on arXiv: Decomposition of mean-field Gibbs distributions into product measures

I’m happy to say that my advisor Ronen Eldan and I recently uploaded a paper to the arXiv under the title “Decomposition of mean-fields Gibbs distributions into product measures” (https://arxiv.org/abs/1708.05859). This is a sister paper of the previous one about exponential random graphs: this one presents a “general framework” and only briefly touches on how it can be used, while the latter is an “applications” paper which uses the framework as a black box to show some neat results.

In this post, I’d like to explain what Gibbs distributions are, what we found out about them, and why all of this is interesting.

This will be an exposition, so I’ll try not to show all the technical details and messy calculations, but some will inexorably wriggle their way in. This shouldn’t stop you from reading the rest of the post, though 🙂 .

Our work is heavily based on my advisor’s prior work concerning Gibb’s distributions, titled “Gaussian-width gradient complexity, reverse log-Sobolev inequalities and nonlinear large deviations” (https://arxiv.org/abs/1612.04346). This post can also be used to familiarize yourself with the background of that paper as well.

Gibbs Distributions

Our paper deals with random vectors $X$ with $n$ elements, where each entry is either $-1$ or $1$. Random vectors with $\{-1,1\}$ entries are very common in the (mathematical) world. Here are three examples of distributions of such vectors.

Coin flips

The vector $X$ can represent $n$ independent flips of a coin, where an entry of $1$ represents a “heads”, and $-1$ represents a “tail”. Suppose that the probability that the coin lands on “heads” is some number $0 \leq p \leq 1$. Then the probability of obtaining some particular vector $x$ that has $H(x)$ heads (and therefore $n-H(x)$ tails) in $n$ coin flips is

$\text{Pr} [X = x] = p^{H(x)}(1-p)^{n-H(x)}.$

If we wish, we can also write this probability as a slightly-cumbersome exponent, using the fact that $p = e^{\log p}$:

$\text{Pr} [X = x] = e^{H(x)\cdot \log p + (n-H(x))\cdot \log (1-p)}.$

In technical terms, this is just a “$\pm 1$ iid Bernoulli random vector with success probability $p$”, which is pretty much the simplest random vector there is.

We can complicate this model a bit by saying that each entry $X_i$ represents the outcome of a different coin, and that coin $i$ has probability $p_i$ to come out heads. In this case, the probability of obtaining a some particular vector $x$ is a product of the form

$\text{Pr} [X = x] = \prod_{i=1}^{n} \begin{cases} p_i & x_i = 1 \\ 1-p_i & x_i = -1 \end{cases}.$

We call this particular distribution a product measure, and say that it is governed by a probability vector $\vec{p}$. These (relatively simple) distributions will play a key role in the things to come.

Exponential random graphs

The exponential random graph model is a way of assigning probabilities to graphs.
I have written about it in greater detail in a previous blogpost, but will give a quick recap now.

In this model, denoted $G_{N}^{f}$, we let $f$ be some function on graphs with $N$ vertices. For example, for a given graph $G$, $f(G)$ may count the number of triangles in $G$; it may return the size of the largest clique in $G$; or it may return 1 if $G$ is a tree and 17 otherwise. Using $f$, we define the following probability distribution on graphs with $N$ vertices:

$\text{Pr} [G_{N}^{f} = G ] = e^{f(G)}/Z,$

where $Z$ is a normalizing constant that is just there to make sure that the probabilities sum up to $1$; it may be written as $Z = \sum_{G} e^{f(G)}.$

Now, a graph with $N$ vertices can have up to ${N \choose 2} = (N^2 - N)/2$ possible edges: One for each pair of vertices. We can then represent a graph by a vector with $n = {N \choose 2}$ entries which take values in $\{-1,1\}$. Each entry represents a possible edge between two vertices in the graph: An entry with the value $1$ indicates that the corresponding edge appears in the graph, while the value $-1$ indicates that the edge does not appear. Thus the probability distribution of the exponential random graph model directly translates into an a probability distribution on vectors.

The Ising model

The Ising model is a model in statistical physics which attempts to capture some of the interesting properties of a magnet. For your leisure and convenience, here is a grossly oversimplified description of the physics of material magnetism; you can skip the following paragraph if all you want is the math.

A solid metal is composed of a bunch of atoms. Each of these atoms has a tiny magnetic moment, which can point either up or down. If all the atoms point upwards together, their magnetic moments add up, and the material is magnetic. The same is true if all the magnetic moments point downwards together. However, if about half of the atoms point upwards and half point downwards, they effectively cancel each other out, and the material is non-magnetic. If this were all there was to it, material-science would be easy; however, the atoms also interact with each other, and this adds a lot of complication to the system. In general, each atom feels the effect of the magnetic moments of its neighbours, and the energy of the system is lower if neighbouring magnetic moments are aligned. Real-life systems tend to be in states with low energy, so a material would “like” to have its all its magnetic moments aligned. But initial conditions + the effects of having a finite temperature usually prevent this.

The Ising model captures the above behaviour in a very simple and elegant way. Our material is modelled as a graph $G$ with $n$ vertices, where each vertex represents an atom. A common example is a (finite) 2D or 3D lattice. In the 2D case, for example, every atom has 4 neighbours (except at the boundary):

The magnetic moments of the atoms are represented by numbers: each vertex gets either a value of $1$ (if the magnetic moment points up) or $-1$ (if the magnetic moments points down). The energy of the system $x$, or the “Hamiltonian of the system” is then defined as

$E(x) = - \sum_{i \sim j \in G} \beta_{ij} x_i x_j,$

where $\beta_{ij}$ is some real number which moderates the strength of interaction between the sites; it is usually controlled by the temperature of the system. If two neighbouring magnetic moments $i$ and $j$ are of equal sign (either both are up or both are down), then the product $x_i x_j$ is $1$ and the energy is reduced (notice the minus sign at the front). But if they are of opposite sign, then $x_i x_j = -1$, and the energy of the system is increased.

Following the rules of statistical physics, the probability of finding a system with low energy should be much higher than the probability of finding a system with high energy. Specifically, if $X$ is an Ising model, the probability of obtaining a particular configuration of the system is

$\text{Pr} [X = x] = e^{-E(x)}/Z,$

where $Z$ is some normalizing constant. In words, this formula means that configurations with low energy are exponentially more likely to appear than configurations with high energy.

This translates directly into a probability distribution of vectors $X \in \{-1,1\}^n$, where $n$ is the number of vertices in the graph $G$: each entry of $X$ corresponds to the direction of the magnetic moment of a vertex in $G$. (Take care not to confuse this with the case of exponential random graphs, where $n$ was the number of possible edges in a graph, and each entry of $X$ corresponded to whether a particular edge was in $G$ or not).

Gibbs distributions

All of the above probability distributions were actually examples of Gibbs distributions. We say that a random vector $X$ follows a Gibbs distributions with Hamiltonian $f$ if the probability for $X$ to be equal to $x$ is just

$\text{Pr} [X = x] = e^{f(x)}/Z,$

where $Z$ is a normalizing constant that is just there to make sure that the probabilities sum up to $1$; it may be written as $Z = \sum_{Z} e^{f(Z)}.$ We denote such a system by $X_{n}^{f}$.

For the case of independent $p$-biased coin flips, the Hamiltonian $f$ was
$f(x) = H(x) \log (p / (1-p)).$ For the case of exponential random graphs, the Hamiltonian $f$ is some graph statistic (e.g the number of triangles in the graph). For the case of the Ising model, $f(x) = \sum_{i \sim j \in G} \beta_{ij} x_i x_j$.

Technical notes (don’t worry if these don’t bother you) :
1) In most sources, Gibbs distributions are defined to have probability proportional to $e^{-f(x)}$ rather than $e^{f(x)}$; there is a minus sign in the exponent. This is purely a matter of convention, and we choose to omit this minus sign in our work.
2) We will mainly want to to talk about asymptotics in this post. For example, we might want to ask how the Ising model on a grid with $n = N \times N$ atoms behaves as the grid size $N$ goes to infinity. Technically, for every $n$ there should be a different function $f$, since $f$ was defined over vectors in $\{-1,1\}^n$. But in this post we will ignore this and think of $f$ as a single function over all possible $n$‘s. This lets us relate $f$ to the size of $n$, so we can say things like “$f(x) < n$ whenever $x$ has $n$ elements”.

Admittedly, the definition by itself doesn’t tell you much about the behaviour of Gibbs distributions, since $f$ was arbitrary, and if we allow $f$ to be arbitrarily complicated, we can get any distribution we desire. That is, for every probability distribution $\pi$ on vectors with $n$-entries, we can always find an $f$ so that $X_{n}^{f} \sim \pi$; simply take

$f(x) = \log{\pi(x)}.$

But that’s ok! Because although an arbitrary distribution can always be obtained by an arbitrary Hamiltonian, we are almost always not interested arbitrary distributions. Rather, the Gibbs distributions that find their way into other applications are usually either simple to state, or simple to calculate, or are considered simple by some other measure. Two of the Hamiltonians given above: the graph triangle-counting function (which is often used in social network theory) and the Ising model (which occupies a large mental space of all statistical physicists) are good examples of this.

Indeed, our main contribution in our paper is to show that there is a notion of complexity, called “gradient complexity”, so that if the Hamiltonian $f$ is “simple” under this notion, then the Gibbs distribution $X_{n}^{f}$ defined by it becomes much more tractable. In what sense is it more tractable? In the sense that it can be well approximated by a mixture of product measures, where the probability vector of each product measure obeys some predefined equation (we’ll talk about what we mean by “mixture” soon). This is in fact the source of the paper’s title, “Decomposition of mean-field Gibbs distributions into product measures”. The adjective “mean-field”, which is a term originating from statistical mechanics, is essentially another name for “low gradient complexity”.

Now you may ask, Why is this a good thing? Why would we even want to decompose Gibbs distributions into mixtures of product measures which obey some predefined equation? The answer is that, in general, Gibbs distributions have so far been quite resilient to analysis, and we don’t understand them very well. You can read the post on exponential random graphs to see an example, but to be brief:

• Although we would very much like to do so, we have a hard time calculating the normalizing constant $Z$.
• Given a Hamiltonian $f$, we have a hard time creating a random sample from $X_{n}^{f}$.
• Given a black box that can generate samples from $X_{n}^{f}$ for some unknown Hamiltonian $f$, we have a hard time figuring out what the original $f$ was.

On the other hand, product measures are much easier to deal with, especially if we know something about their probability vectors. Of course, as stated above, we don’t claim that the distribution is close to any particular product measure, but rather to a “mixture of product measures”. Suppose that $\rho$ is a probability measure on vectors in $[-1,1]^n$. A distribution is called a “$\rho$-mixture”, denoted $X(\rho)$, if it can be randomly sampled as follows: first, generate a random vector $z$ according to the measure $\rho$. Then, generate a random vector according to the product measure whose probability vector has expectation $z$. In other words, we first generate random coin-biases using $\rho$ (by randomly generating the coins’ expectations), and then we flip $n$ independent coins, one for each generated bias.

Advanced note: As usual, unless we know something about $\rho$, this definition doesn’t give us a lot of information: every distribution can be expressed as a $\rho$-mixture, for $\rho$ which mimics the original distribution by simply has all its mass on $\{0,1\}^n$ singletons. Our results are meaningful because they show that $\rho$ has a very specific structure called “small tilt decomposition” which gives it useful properties, but we won’t go into that in this post.

Our main result can thus be very, very, very roughly written as

Theorem 1, rough statement: If the Hamiltonian $f$ is nice enough (i.e has low gradient complexity + some other smoothness conditions) then there exists a measure $\rho$ so that:

1. $X^f_n$ can be approximated by a $\rho$-mixture.
2. Most of the mass of $\rho$ rests on vectors $x$ which satisfy

$\frac{1}{n} \| x - \tanh(\nabla f(x)) \|_1 \to 0$

as $n \to \infty$, where the $\tanh$ is applied entrywise to the entries of $\nabla f$. Here $\|x\|_1$ is the one-norm of the vector $x$, and stands for

$\|x\|_1 = \sum_{i=1}^{n} |x_i|.$

We have already talked a bit about the first item in the above theorem, so let’s inspect the second one. I haven’t told you yet what $\nabla f(x)$ means in this context, but for now it’s enough to think of it as an analogue to the real-differentiable gradient of a function: it is a vector with $n$ elements, where each coordinate is some simple function of $f(x)$. The important thing to see is that this equation tells us something about how the probability vectors of the individual product measures behave: for almost all the vectors which will appear under $\rho$, we know that they are close to satisfying the fixed point equation

$x = \tanh(\nabla f(x)).$

If we can solve this equation, then we can say something meaningful about how $x$ looks like, and hence how the product measure looks like.

I say “are close to satisfying the fixed point equation…”, because there can still be some error: all we require is that $\| x - \tanh(\nabla f(x)) \|_1$ is asymptotically smaller than $n$, something which is often written as

$\| x - \tanh(\nabla f(x)) \|_1 = o(n).$

This means that almost every entry of $x$ is close to the corresponding entry of $\tanh(\nabla f(x))$; indeed, if the distance between them as vectors is $o(n)$, then almost every individual entry is $o(1)$. Of course, this may not always be an easy equation to solve.

As we will soon see, there is a technical and precise way to define what “nice” means in this context. But before we do that, I’d first like to show some intuition as to what sort of Hamiltonians might be decomposable into product measures, and what sort of Hamiltonians might not.

Intuitively, we might think that since product measures have independent entries, mixtures of product measures also still contain a lot of independence. That is, if I know that a certain element in my vector is, say, 1, this shouldn’t affect the probability of other vertices being 1 by too much. Thus, Hamiltonians for which there is not too strong a dependence between $X$‘s entries should be decomposable, while Hamiltonians with very strong dependence between entries should not. Of course, for a Gibbs distribution to be interesting it must have at least some dependence between the entries of $X$, otherwise it is just a product measure.

An example for a Hamiltonian which does not satisfy this requirement is the 2D, 4-neighbor Ising model described above. For suppose that I know that the magnetic moment of a certain atom is pointing upwards. This has a great impact on its four neighbors: I am now certain about 1/4th of the forces acting on the neighboring atoms, and can be much more confident that they will point upwards as well. Thus there is a (relatively) strong dependence between adjacent nodes in this model (further, this dependence propagates: since we know that the neighbors are more likely to point upwards, their neighbors are also more likely to point upwards, and so on; but the effect is felt strongest in the immediate vicinity of the initial site).

On the other hand, consider a similar 2D Ising model where each atom is now connected to a large neighborhood: all pairs of sites are considered neighbors if they are within distance $R$ from each other, where $R$ is not a constant but increases as the size of the lattice grows.

In this model, the effect of any individual atom on its neighbors is minuscule: every magnetic moment is influenced by a large number of other atoms, and the dependence of any atom on any specific individual neighbor decreases to 0 as $n \to \infty$. We can therefore expect neighboring sites to be “almost independent”, something which hopefully can be captured by a mixture model of product measures.

Another way to see the difference between the models is as follows. In the regular 2D model, where each site only has 4 neighbors, I can pick a site, and fix its 4 neighbors to be either all up, or all down. If I fix them to point upwards, then with relatively high probability, the atom in the middle will also point upwards; likewise if I fix them to point downwards. The whole affair only deals with 5 atoms, so very roughly speaking, I can fix about $n / 5$ of the sites to point in any arbitrary directions I wish (there are problems at the interface between such clumps, but let’s not worry about that, this is just for intuition).

On the other hand, if each site is influenced by many other sites, as in the second model we described, things are much “smoother”. If we wish to control the direction of a particular site, we have to fix a lot of magnetic moments in the ball around that site. Since the radius of the ball increases with $n$, this means that we can only arbitrarily control $o(n)$ sites.

Thus, for the 4-neighbor model, we can create probable vectors where about $n / 5$ of the sites point in arbitrary directions, but for the large-ball model, only $o(n)$ of them can be made arbitrary. In this sense, the latter model is much smoother than the former. Our definition of “Gaussian Width” will utilize this observation.

Anyway, back to the theorem statement; to fully understand it, I still owe you several things:

1. What does it mean for $f$ to be “nice”?
2. What does it mean to be approximated by a $\rho$-mixture?
3. What does $\nabla f(x)$ mean for a function $f$ that is normally defined on $\{-1,1\}^n$?

I’ll now give answers to these three questions. The details will be a bit technical (but only a bit!). If you are fine with the above theorem formulation as it is, though, you can skip on ahead to the next section of this post (titled “An application”).

1) We start with what it means for $f$ to be nice, and for that we need define what it means to “differentiate” a function whose domain is $\{-1,1\}^n$. For a coordinate $i$, define the $i$-th partial derivative of $f$ as

$\partial_i f(y) = \frac{1}{2}(f(y_1,\ldots, y_{i-1}, 1, y_{i+1}, \ldots,y_n) - f(y_1,\ldots, y_{i-1}, -1, y_{i+1}, \ldots,y_n)).$

In other words, we look at the change in the value of $f$ when we flip just a single coordinate. This lets us define two important quantities. The first is the discrete gradient of $f$:

$\nabla f(y) = (\partial_1 f(y), \ldots, \partial_n f(y)),$

which is just a straightforward analog of the real-differentiable gradient. The second is the Lipschitz constant of $f$:

$\text{Lip}(f) = \max_{i\in[n],y\in\{-1,1\}}|\partial_i f(y)|.$

This is reminiscent of the good ol’ Lipschitz constant we are used to having in standard analysis, since by the triangle inequality, for every $x,y \in \{-1,1\}^n$ we have

$|f(x) - f(y)| \leq \text{Lip}(f) \|x-y\|_1.$

The Gaussian Width of a set $K \subseteq \mathbb{R}^n$ is defined as

$\mathbf{GW}(K) = \mathbb{E}[\sup_{x\in K}\langle x, \Gamma \rangle],$

where $\Gamma \sim N(0,\text{Id})$ is a standard Gaussian vector in $\mathbb{R}^n$.

In other words, we do the following. First, we generate a random Gaussian vector $\Gamma$, that is, a vector where each of its entries is a $N(0,1)$ standard good-ol’-Gaussian random variable. Then we look inside the set $K$ for the point with the highest correlation with $\Gamma$. We then take an average of this over the distribution of $\Gamma$.

The name “Gaussian Width” comes from the fact that this quantity, in a way, represents the “width” of the set $K$, only normalized by the factor of a Gaussian. The gradient complexity of a function is then just defined as

$\mathcal{D}(f) = \mathbf{GW}(\{\nabla f(y) : y \in \{-1,1\}^n \} \cup {0}).$

In words, it is the Gaussian Width of the range of the discrete gradient of $f$.

So $\mathcal{D}(f)$ is the average highest correlation of a random Gaussian with the possible gradients $\nabla f$. We can relate this to the intuition we described earlier, when we contrasted the non-smoothness of the 4-neighbor Ising model with the smoothness of the large-ball Ising model.

Suppose you already generated the vector $\Gamma$. For the 4-neighbor model, look at the largest $n / 5$ entries in absolute value of $\Gamma$. About half of them will be positive, and about half of them will be negative. But because we can arbitrarily control $n / 5$ of the entries of $\nabla f$, there exists a state which highly correlates with those $n / 5$ entries of $\Gamma$: just make sure that the signs align. Since the entries of $\Gamma$ are independent, the rest of the entries won’t contribute much either positively or negatively to the correlation, and we get that $\mathcal{D}(f)$ is of order of magnitude $n$.

On the other hand, this trick cannot be repeated with the large-ball model, since there we can only arbitrarily fix $o(n)$ entries. While this doesn’t by itself mean that the large-ball model has low complexity – this still needs to be proven formally – it gives you a hunch of why it doesn’t have high complexity: the gradient is too “smooth” to have a large correlation with random noise (which is just what a Gaussian vector is).

So what is a “nice” function? It is a function for which:

1. $\mathcal{D}(f)$ is asymptotically smaller than $n$, i.e $\mathcal{D}(f) / n \to 0$ as $n \to \infty$. Intuitively, this means that the set of possible gradients of $f$ doesn’t take up too large a volume.
2. $f$ has a small Lipschitz constant; e.g $\text{Lip}(f)$ does not depend on $n$. This means that $f$ doesn’t change by too much as you change its input.
3. $\nabla f$ has a small Lipschitz constant (in the metric-space sense); e.g, one that does not depend on $n$. This means that $\nabla f$‘s norm doesn’t change by too much as you change its input.

As you can see, these conditions are all different ways of saying that $f$ basically can’t change too much, or too quickly. In this sense, it is “well behaved”.

2) When we say that $X_n^f$ is approximated by a $\rho$-mixture, we mean that we can generate vectors both under the law of $X_n^f$ and under the law of the $\rho$-mixture in such a way, so that in expectation, the distance between them is small. The formal definition goes like this: Let $\delta > 0$. A random variable $X$ is called a $(\rho, \delta)$-mixture if there exists a coupling between the $\rho$-mixture $X(\rho)$ and the original random vector $X$ so that

$\mathbb{E}\|X(\rho) - X\|_1 \leq \delta n.$

Being well approximated corresponds to having a $(\rho, \delta)$-mixture with $\delta \to 0$ as a function of $n$; in other words, the expected distance between vectors of the two distributions is asymptotically smaller than $n$. This implies that almost all of the time, the differences between almost all of the individual coordinates tend to $0$.

3) If you are especially keen eyed, you may have noticed a discrepancy between the equation $\| x - \tanh(\nabla f(x)) \|_1 = o(n)$ and the definition of the discrete gradient $\nabla f$: In the equation, the entries of $x$ take values in the continuous cube $[-1,1]^n$, while the discrete gradient $\nabla f$ is only defined over the discrete cube $\{-1,1\}^n$. How can we calculate $\nabla f$ for continuous vectors? It turns out that there are two natural ways to extend a discrete function to the entire continuous cube; I’ll describe them here only very briefly without going into the details, so don’t worry if you don’t know all of these terms. The first is to take the discrete Fourier decomposition of $f$ into monomials, and just plug in fractional values into the resulting polynomial formula. The second is to treat the discrete values of $f$ as boundary conditions for solving Laplace’s partial differential equation, and progressively solve the problem in higher and higher dimensions. As it happens, the two methods are actually equivalent! So you can use your favorite one to extend $\nabla f$.

We can now finally state Theorem 1 in a more precise version.

Theorem 1, in almost all its glory: Let $f$ be a Hamiltonian and denote

$D = \mathcal{D}(f)$

$L_1 = \max \{1, \text{Lip}(f)\}$

$L_2 = \max \{1, \max_{x\neq y} \dfrac{\| \nabla f(x) - \nabla f(y)\|_1}{\|x-y\|_1} \}.$

Denote by $\mathcal{X}_f$ the set

$\mathcal{X}_f = \{x \in \{-1,1\}^n: \| x - \tanh(\nabla f(x)) \|_1 \leq 600 e^{2L_1} L_2^{3/4} D^{1/4} n^{3/4} \}.$

Then there exists a measure $\rho$ such that $X_n^f$ is a $(\rho, 20\frac{D^{1/4}}{n^{1/4}})$-mixture, with

$\rho (\mathcal{X}_f) \geq 1 - 20\frac{D^{1/4}}{n^{1/4}}.$

Here the set $\mathcal{X}_f$ is the set of “good vectors”, which are close to obeying the fixed point equation $x = \tanh(\nabla f(x))$. If you recall, we said that “nice” functions are those for which $D = o(n)$, and $L_1$ and $L_2$ are constant. In this case, the error term in the definition of $\mathcal{X}_f$ is $o(n)$, the $\delta$ term in definition of the $(\rho,\delta)$-mixture goes to $0$, and the mass of the good set $\mathcal{X}_f$ goes to $1$ as $n \to \infty$. The formulation above also shows that $L_1$ and $L_2$ don’t have to be constant, but are allowed to grow somewhat with $n$ while still retaining a good bound.

The above formulation is only in “almost all its glory” because, we actually show that $X_n^f$ is something stronger than a $\rho$-mixture – we show that it is a tilt decomposition – but for that you are going to have to read the paper itself :).

An application

Now that we have a hammer, it’s to time start looking for nails. The most straightforward way to apply our theorem is to take your favorite Gibbs distribution and start calculating the gradient complexity and Lipschitz constants for its Hamiltonian. If these are small enough, hurray! You can try to solve the equation $\|x - \tanh (\nabla f(x))\|_1 = o(n)$ in order to say something interesting about the product measures in the decomposition.

Let’s see an example of this: The Curie-Weiss Ising model. In this model, we have a collection of $n$ atoms which all interact with each other; that is, every magnetic moment feels the effect of every other magnetic moment. The underlying graph for this model is then the complete graph. Admittedly, the model doesn’t capture physical reality very well: it is not really possible, in large blocks of matter, for all the atoms to interact with one another. However, it is useful to study such simple things, in order to gain insights into the behaviour of more complicated models.

The Hamiltonian $f$ is written as

$f(x) = \beta \frac {1}{n} \sum_{i\neq j} x_i x_j.$

Note that there is a $1/n$ term in front of the sum; the rationale for this is that the force on a single atom cannot explode to infinity, so in fact each atom only feels the average of all its neighbors. The term $\beta$ is often called the “inverse temperature”; if we were to attach a physical meaning to it, it would be proportional to $1/T$, where $T$ is the temperature. A temperature tending to $0$ will give a gigantic $\beta$ which puts a strong emphasis on the interaction between atoms, while a temperature tending to $\infty$ will means that the interaction between atoms becomes negligible.

So, the first task is check whether the Curie-Weiss Hamiltonian is “nice”. This is not very hard to do once you treat the problem formally: it turns out that you can express the Hamiltonians of Ising models using matrices: any Hamiltonian of the form

$f(x) = \sum_{i\neq j} a_{ij} x_i x_j$

can be written as

$f(x) = x^T A x$

where $A_{ij} = (a_{ij})$ is a symmetric matrix called the “interaction matrix”. Under this notation, the gradient of $f$ can be calculated to be

$\nabla f(x) = Ax,$

so all we have to do is check the Gaussian Width and Lipschitz constants of a linear operator. For the case of the Curie-Weiss model, this operator is luckily very simple: all the interaction coefficients $a_{ij}$ are the same and are equal to $\beta/n$. The interaction matrix $A$ is then just

$A = \begin{bmatrix} 0 & \beta/n & \beta/n & \ldots & \beta/n \\ \beta/n & 0 & \beta/n & \ldots & \beta/n \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \beta/n & \beta/n &\beta/n & \ldots & 0 \end{bmatrix}.$

We won’t do the calculations here, but analyzing this matrix shows that the Curie-Weiss Hamiltonian indeed satisfies our conditions for being “nice”. So we can apply our theorem, and try so see how the individual product measures of the decomposition behave.

Now, as it turns out, in this particular case, instead of grinding our gears and trying to solve the approximate fixed point equation $\|x - \tanh (Ax)\|_1 = o(n)$, we can opt to solve the exact fixed point equation

$x = \tanh (Ax).$

The solutions to the exact equation will tell us basically everything we need to know about the approximate ones. Further, it can be shown that not much is lost if we assume that the matrix $A$ has $\beta/n$ on the diagonal as well, instead of $0$‘s. In this case, the above equation reduces to a scalar equation, involving no matrices at all! This is because then every entry $x_i$ of a solution to the vector equation $x = \tanh(Ax)$ satisfies

$x_i = \tanh(\sum_{j=1}^{n}\frac{\beta}{n}x_j),$

and so all the $x_i$‘s are equal to each other! Every solution of this equation is then of the form $x = (a,a,\ldots,a)$ for some real number $a$. We then only need to solve the real-valued equation

$a = \tanh(\beta a).$

As the following image shows, for $\beta < 1$, there is only one solution, which is $a = 0$.

Let’s interpret this result. The vectors $x$ which solve the fixed point equation represent the expectation of a product measure. If every coordinate of the expectation is 0, then every coordinate of the product measure has a probability $1/2$ of getting the value $1$, and a probability $1/2$ of getting the value $-1$. That is, there is absolutely no preference between magnetic moments which point up, and magnetic moments which point down. Whenever we sample a system, we can expect about half the sites to point up, and about half the sites to point down.

This is what is called in physics an “unordered state”, which, magnetism-wise, corresponds to “no magnetism”. Recalling that $\beta$ represents the inverse temperature, a small value of $\beta$ corresponds to a high temperature, so this can also be stated as follows: when the temperature is high, there can be no magnetism. This appeals to the physical intuition that when it is very very hot, the thermal fluctuations which randomly flip the magnetic moments are so large that they do not allow the system to settle down in any ordered state.

On the other hand, when $\beta > 1$, there can be more than one solution: in fact, two new solutions appear, one for positive $a$, and one for negative $a$:

These solutions tend towards $1$ and $-1$ as $\beta \to \infty$. What do these solutions mean? A value of $a = 1/2$, for example, indicates that for each site, there is a probability of $3/4$ for the atom to point upwards, and a probability of $1/4$ for the atom to point downwards; that way, the expectation is indeed $1/2$. In this case, about $3/4$-ths of the magnetic moments will point upwards, and we have an ordered state: a magnet. This happens at low temperatures, and appeals to the physical intuition that when it is very cold, there are almost no thermal fluctuations, and the system can settle down in a state of low energy – that is, one where most of the magnetic moments point in the same direction.

Final note

Our paper contains all formalities and proofs and explanations to everything that was left unexplained above (or at least, references to proofs and explanations). I hope that this post sparked your interest enough so that you would take a look at the details. The application of our main theorem to the Curie-Weiss Ising model is only a toy example; perhaps there are many other applications which are waiting to be picked apart. For instance, we ourselves used this framework to say something interesting about exponential random graphs. But that’s already a story for another post…