<figure><img src="image_1.png" /></figure>

HAL
open science

VARIATIONAL INFERENCE OF A TWO-STATE
SYSTEM

Weston D Viles

► To cite this version:

Weston D Viles. VARIATIONAL INFERENCE OF A TWO-STATE SYSTEM. 2025. hal-05231514

HAL Id: hal-05231514

https://hal.science/hal-05231514v1

Preprint submitted on 30 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# VARIATIONAL INFERENCE OF A TWO-STATE SYSTEM

WESTON D. VILES

**ABSTRACT.** A probability distribution over the unit hypercube encodes incomplete information on the configuration of a system of interacting entities that exist in one of two states represented by a vector indicator random variables. We infer the probability distribution of the sum of jointly distributed indicator random variables under arbitrary statistical dependence from the first principles of information and probability theory. We provide a tractable, low-rank approximation to the exact probability distribution with a corresponding expression for the approximation error. Our method is applied to the partition function of the ferromagnetic Ising model to derive a deterministic approximation thereof as a demonstration of its utility.

**Key phrases:** sum of dependent indicator random variables, principle of maximum entropy, variational inference, ferromagnetic Ising model

## 1. BACKGROUND

Any two-state system of $n \in \mathbb{Z}_+ = \{1, 2, \dots\}$ interacting entities in equilibrium may be represented by a member of the set $Q = \{Q : \{0, 1\}^n \to \mathbb{R}_+ \cup \{0\} \mid \|Q\|_1 = 1\}$ of probability distributions over the unit hypercube, where $\mathbb{R}_+ = \{x \in \mathbb{R} \mid x > 0\}$. Let the support $\Omega = \{\mathbf{x} \in \{0, 1\}^n \mid Q(\mathbf{x}) > 0\}$ of $Q$ be partitioned by the sets $\Omega_x = \{\mathbf{x} \in \{0, 1\}^n \mid \|\mathbf{x}\|_1 = x\}$ for $x \in \mathbb{N}_n = \{0, 1, \dots, n\}$. The induced probability distribution

$$P(x) = \sum_{\mathbf{x} \in \Omega_x} Q(\mathbf{x})$$

is that of the sum $X = X(\mathbf{x}) = x_1 + \dots + x_n$ of the components of $\mathbf{x} = (x_1 \dots x_n)' \in \{0, 1\}^n$. The prior distribution $Q'(\mathbf{x}) \propto 1$ of $Q$ for $\mathbf{x} \in \Omega$ induces the prior distribution $P'(x) \propto |\Omega_x|$ of $P$ for $x \in \Upsilon = \{x \in \mathbb{N}_n \mid P(x) > 0\}$ in the support of $P$. If the system is ergodic then $P'(x) \propto \binom{n}{x}$ for $x \in \mathbb{N}_n$.

The probability distribution of a sum of jointly distributed indicator random variables, if not a binomial distribution, is classically approximated by a binomial or Poisson distribution with error quantified in variation distance: the quintessential metric in Kolmogorov's [23] set and measure theory of probability. The error associated with an appropriately parameterized binomial or Poisson approximation may be governed by imposing strict conditions on the joint distribution of the indicator random variables, e.g., weak dependence or positive/negative relations, etc. [1, 2] The structure and rigidity of these one-parameter families, as well as that of their fixed-parameter extensions, hinders their utility in scientific applications. See [27, 9, 5, 3, 35] and references therein for original results and extensions thereof.

<footnote>Date: August 30, 2025.</footnote>
The tenet that probability distributions encode incomplete information about a
system underlies the fundamental work of Shannon, [39] Jaynes, [16, 17] and Kull-
back. [26, 24] *Information entropy* quantifies the incomplete information encoded
by a probability distribution over states of the system as the average information
acquired from a single observation of the system. The *principle of maximum en-
ropy* posits the statistical representation of a system as that which maximizes
information entropy subject to macroscopic constraints on the system. Jaynes'
theory of probability [18] formalizes inference with incomplete information in con-
formity with the scientific method. Further development of the theory delivered the
*principle of maximum relative entropy* [13] and, ultimately, the *exponential family*
of probability distributions of mathematical statistics. [8]

Variational inference is the method of implicitly defining a quantity of interest as the solution to an optimization problem. [44] We infer the probability distribution P with the variational approach prescribed by the principle of maximum entropy and provide a low-rank approximation to the exact distribution of arbitrary precision. The utility of our methodology is demonstrated with the derivation of a deterministic approximation to the partition function of the ferromagnetic Ising model [28, 15] on general simple graphs.

2. ON THE SUM OF INDICATOR RANDOM VARIABLES

The support $\Upsilon \subseteq \mathbb{N}_n$ of the probability distribution $P$ corresponds to the support $\Omega \subseteq \{0, 1\}^n$ of $Q$ since $P(x) > 0$ for some $x \in \mathbb{N}_n$ if and only $Q(x) > 0$ for some $x \in \Omega_x$. We suppose that $0 \in \Upsilon$, viz., $0 \in \Omega$, without loss of generality and define $\Upsilon_0 = \Upsilon \setminus \{0\}$ accordingly.

Let $\mathfrak{e}_k : \{0, 1\}^n \to \mathbb{N} = \{0, 1, \dots\}$ be the $k$th-order elementary symmetric function for $k \in \mathbb{N}$. Since $\mathbf{x} \in \{0, 1\}^n$ then $\mathfrak{e}_k(\mathbf{x}) = \binom{X}{k}$ and

$$
\begin{align*}
\sum_{x \in \Upsilon} \binom{x}{k} P(x) &= \sum_{x \in \Upsilon} \binom{x}{k} \sum_{\mathbf{x} \in \Omega_x} Q(\mathbf{x}) \\
&= \sum_{x \in \Upsilon} \sum_{\mathbf{x} \in \Omega_x} \mathfrak{e}_k(\mathbf{x}) Q(\mathbf{x}) = \sum_{\mathbf{x} \in \Omega} \mathfrak{e}_k(\mathbf{x}) Q(\mathbf{x}),
\end{align*}
$$

that is, $E_P(\binom{X}{k}) = E_Q \mathfrak{e}_k(\mathbf{x})$ for $k \in \mathbb{N}_n$, where $E_P$ and $E_Q$ are the expectation operators with respect to $P$ and $Q$, respectively. The population moment

$$
E_P \binom{X}{k} = \sum_{\substack{J \in \mathbb{N}_n \\ |J| = k}} E_Q \left( \prod_{j \in J} x_j \right)
$$

with respect to $P$ is thus a sum of corresponding low-order population moments
with respect to $Q$ for $k \in \mathbb{N}_n$. The set of $|\Upsilon| \in \mathbb{N}_n$ moments $E_P(\binom{X}{k}) = E_Q \mathfrak{e}_k(\mathbf{x})$ are
linearly independent for $k \in \Upsilon$, i.e., if $\Upsilon \subset \mathbb{N}_n$ and $k' \in \mathbb{N}_n \setminus \Upsilon$ then

$$
E_P \binom{X}{k'} = \sum_{\substack{x \in \Upsilon \\ x \ge k'}} \binom{x}{k'} P(x) = E_P \binom{X}{\min\{k \in \Upsilon \mid k > k'\}}
$$

if {$k \in \Upsilon \mid k > k'$} $\neq \emptyset$ and $E_P(\binom{X}{k'}) = 0$ otherwise. The support $\Upsilon \subseteq \mathbb{N}_n$ of $P$, if unknown a priori, may thus be determined through the linear (in)dependence of population moments.
Let $\mathcal{H}_\Upsilon = \{h : \Upsilon \mapsto \mathbb{R}_+\}$ be the set of functions with support $\Upsilon \subseteq \mathbb{N}_n$ and consider the constrained optimization problem

$$ \min_{h \in \mathcal{H}_\Upsilon} \text{KL}(h \| P') \\ \text{subject to} \quad E_h\binom{X}{k} = E_p\binom{X}{k} \text{ for } k \in \Upsilon, $$

in which KL ($h \| P'$) is the Kullback-Leibler divergence [26] from the prior distribution $P'$ of $P$ to the probability distribution $h \in \mathcal{H}_\Upsilon$. Let $h \in \mathbb{R}_+^{n+1}$ be zero-indexed and $\theta \in \mathbb{R}^n$ be one-indexed with respective components $h_x \in \mathbb{R}_+$ for $x \in \Upsilon$ and $h_x = 0$ otherwise and $\theta_k \in \mathbb{R}$ for $k \in \Upsilon_0$ and $\theta_k = 0$ otherwise. Let $\theta_0 \in \mathbb{R}$ and define $L(\mathbf{h}, \boldsymbol{\theta}, \theta_0)$ as the Lagrangian

$$ L(\mathbf{h}, \boldsymbol{\theta}, \theta_0) = \sum_{x \in \Upsilon} \left[ h_x \log \left( \frac{P'(x)}{h_x} \right) + \sum_{k \in \Upsilon} \theta_k \binom{x}{k} (h_x - P(x)) \right] $$

corresponding to the constrained optimization problem above.

The components of the stationary point of the Lagrangian $L(\mathbf{h}, \boldsymbol{\theta}, \theta_0)$ corresponding to $\mathbf{h}$ and $\theta_0$ are the solution to $\mathbf{0} = \nabla_{\mathbf{h}} L(\mathbf{h}, \boldsymbol{\theta}, \theta_0)|_{\Upsilon}$ and $0 = \frac{d}{d\theta_0} L(\mathbf{h}, \boldsymbol{\theta}, \theta_0)$, where $\mathbf{v}|_{\Upsilon}$ is the subvector of $\mathbf{v} \in \mathbb{R}^n$ restricted to indices in $\Upsilon \subseteq \mathbb{N}_n$. The solution $h_\theta : \Upsilon \mapsto \mathbb{R}_+$ takes the functional form

$$ h_{\theta}(x) = \begin{cases} |\Omega_x| \exp\{\sum_{k \in \Upsilon_0} \theta_k \binom{x}{k} - \Psi_h(\theta)\} & \text{for } x \in \Upsilon \\ 0 & \text{otherwise,} \end{cases} $$

where $\Psi_h(\theta) = \log Z_h(\theta)$ is the cumulant function, i.e., log partition function, [44] of the probability distribution $h_\theta$ and

$$ Z_h(\theta) = \sum_{x \in \Upsilon} |\Omega_x| \exp \left\{ \sum_{k \in \Upsilon_0} \theta_k \binom{x}{k} \right\} $$

is the partition function of $h_\theta$. The functions $t_k : \mathbb{N}_n \mapsto \mathbb{N}$ defined as $t_k(x) = \binom{x}{k}$ for $k \in \mathbb{N}_n$ are sufficient statistics of the probability distribution $h_\theta$ that belongs to the exponential family of probability distributions.

The reduced Lagrangian $L'(\theta) = L(h_\theta, \theta, 1 - \Psi_h(\theta))$ takes the expression

$$ \begin{aligned} L'(\theta) &= \sum_{x \in \Upsilon} \left[ h_{\theta}(x) \log \left( \frac{|\Omega_x|/|\Omega|}{h_{\theta}(x)} \right) + \sum_{k \in \Upsilon_0} \theta_k \binom{x}{k} (h_{\theta}(x) - P(x)) \right] \\ &= -\log |\Omega| + \Psi_h(\theta) - \sum_{x \in \Upsilon} P(x) (\log h_{\theta}(x) + \Psi_h(\theta)) \\ &= \text{KL}(P \| h_{\theta}) + H(P) - \log |\Omega| \end{aligned} $$

upon partial identification of the stationary point of $L(\mathbf{h}, \boldsymbol{\theta}, \theta_0)$ over $\mathbf{h}$ and $\theta_0$, where $H(P)$ is the information entropy of $P$. Since the sufficient statistics $t_k$ for $k \in \mathbb{N}_n$ are linearly independent and

$$ E_P t_k = P(k) + \sum_{x \in \Upsilon \setminus \{k\}} \binom{x}{k} P(w) \geq P(k) > 0 $$

for $k \in \Upsilon$ then there exists $\vartheta \in \mathbb{R}^n$ that solves the system $\mathbf{0} = \nabla_{\boldsymbol{\theta}} L'(\boldsymbol{\theta})|_{\Upsilon_0}$ with $\vartheta_k = 0$ for $k \neq \Upsilon_0$.
The Hessian matrix

$$ \nabla_{\theta}^{2} L'(\theta)|_{\Upsilon_{0} \times \Upsilon_{0}} = \mathrm{Cov}_{h_{\theta}}(\mathbf{t})|_{\Upsilon_{0} \times \Upsilon_{0}} $$

is positive definite, where $\mathbf{t} = (t_1 \dots t_n)' : \mathbb{N}_n \to \mathbb{N}^n$ is the vector of sufficient statistics and $M|_{\Upsilon_0 \times \Upsilon_0}$ is the $|\Upsilon_0| \times |\Upsilon_0|$ submatrix of $M \in \mathbb{R}^{n \times n}$ restricted to indices in $\Upsilon_0 \subseteq \mathbb{N}_n$. The stationary point

$$ (1) \qquad \vartheta = \arg \min_{\theta \in \mathcal{R}_{\Upsilon_0}^n} \mathrm{KL}(P\|h_{\theta}) $$

of $L'(\theta)$ is therefore unique, where

$$ \mathcal{R}_{\Upsilon_0}^n = \bigtimes_{k=1}^n \begin{cases} \mathbb{R} & \text{if } k \in \Upsilon_0 \\ \{0\} & \text{otherwise.} \end{cases} $$

The convex function $\mathrm{KL}(P\|h_{\theta})$ may be minimized numerically with a gradient-based approach, although full specification of $h_{\vartheta} = P$ requires the evaluation of $\sum_{k \in \Upsilon_0} \binom{n}{k}$ population moments over $\Omega$ with respect to $Q$ thereby necessitating a low-rank approximation to $h_{\vartheta}$ for tractable application of the method.

Let $v_0 \subseteq \Upsilon_0$ and define

$$ \mathfrak{h}_{\phi}(w) = \begin{cases} |\Omega_x| \exp \{ \sum_{k \in v_0} \phi_k t_k(x) - \psi_{\phi}(\phi) \} & \text{for } x \in \Upsilon \\ 0 & \text{otherwise,} \end{cases} $$

in which $\phi \in \mathcal{R}_{v_0}^n$ and $\psi_{\phi}(\phi)$ is the cumulant function of the probability distribution $\mathfrak{h}_{\phi}$. Define

$$ \varphi = \arg \min_{\phi \in \mathcal{R}_{v_0}^n} \mathrm{KL}(P\|\mathfrak{h}_{\phi}), $$

so that the average additional information per single observation of the system modeled by $\mathfrak{h}_{\varphi}$ in place of $h_{\vartheta}$ increased by

$$ \mathrm{KL}(h_{\vartheta}\|\mathfrak{h}_{\varphi}) = (\varphi - \vartheta)|_{\Upsilon_0} \cdot (E_{\mathfrak{h}_{\varphi}}\mathbf{t} - E_{h_{\vartheta}}\mathbf{t})|_{\Upsilon_0} + O(\|\varphi - \vartheta\|^2), $$

where the positive constant associated with the big-O term is guaranteed by a positive definite Hessian matrix.

### 3. APPLICATION: THE FERROMAGNETIC ISING MODEL

Let $G = (\mathbf{V}, \mathbf{E})$ be a simple graph with $n = |\mathbf{V}|$ vertices and $m = |\mathbf{E}|$ edges. The ferromagnetic Ising model is the probability mass function

$$ (2) \qquad \mathrm{Pr}(\mathbf{z} | G, \alpha, \beta) \propto \exp \left\{ \alpha \sum_{w \in \mathbf{V}} z_w + \beta \sum_{(u,v) \in \mathbf{E}} \beta_{uv} z_u z_v \right\} $$

with support $\mathbf{z} \in \{-1, 1\}^n$ and parameterized by $\alpha \in \mathbb{R}$ and $\beta \in \mathbb{R}_+$. The model is a simplified representation of ferromagnetism that was proposed by Lenz [28] and investigated by Ising [15] for a description of an order-disorder phase transition. See [34, 31, 4] for the history and relevance of the Ising model in the statistical mechanics literature.

As one of the simplest representations of a two-state system of interacting entities suitable for modeling nature, the Ising model has since been of interdisciplinary interest, [30] including in applications of neuroscience [12, 37, 46] and artificial
intelligence. [14, 38] Its utility as a statistical model, however, is impeded by the unyielding analytic and numeric intractability [45, 42, 11] of its partition function

$$Z(G | \alpha, \beta) = \sum_{\mathbf{z} \in \{-1, 1\}^n} \exp \left\{ \alpha \sum_{w \in \mathbf{V}} z_w + \beta \sum_{(u,v) \in \mathbf{E}} z_u z_v \right\}$$

on general graphs. The partition function Z is the normalizing constant implicit in Equation (2) and characterizes the probability distribution by its moment sequence.

Analytic solutions exist on planar graphs, [21, 10] including the one- and two-dimensional lattices [15, 33, 22] and trees [6], but are elusive on general simple graphs. Existing methods of deterministic approximation of Z are restricted to certain graphical regimes, including bounded degree graphs [40, 41] or are otherwise based upon correlation decay and zeros of the partition function in the complex plane. [32, 29, 36] We provide a low-rank approximation to the partition function of the ferromagnetic Ising model on general graphs based upon the methods described in the preceding section: a result that translates to a tractable numerical procedure.

**3.1. Probabilistic formulation.** Let $\mathcal{G} = \{g = (\mathbf{V}, e) \mid e \subseteq \mathbf{E}\}$ be the set of spanning subgraphs of $G$. The partition function Z may be expressed as

$$(3) \qquad Z(G | \alpha, \beta) = (2 \cosh \alpha)^n (2 \cosh \beta)^m \mathcal{Z}(G | \xi, \zeta),$$

where $\xi = \tanh \alpha$ and $\zeta = \tanh \beta$, and the function

$$\mathcal{Z}(G | \xi, \zeta) = \frac{1}{|\mathcal{G}|} \sum_{g \in \mathcal{G}} \xi^{X(g)} \zeta^{Y(g)},$$

in which $X(g)$ and $Y(g)$ are respectively the number of vertices of odd degree and the number of edges in the spanning subgraph $g \in G$. Up to the factor $1/|\mathcal{G}| = 1/2^m$ that causes $\mathcal{Z}$ to be an expected value over a finite population, i.e., the probabilistic formulation, the above expression is proportional to the combinatorial formulation [43, 20] or high-temperature expansion [31, 19] of the partition function in the statistical mechanics literature.

Let $\mu(x, y) = |\mathcal{G}_{xy}|/|\mathcal{G}|$ be the joint distribution of $X = X(g)$ and $Y = Y(g)$, where $\mathcal{G}_{xy} = \{g \in \mathcal{G} \mid X(g) = x, Y(g) = y\}$. The expression for $\mathcal{Z}$ thus reduces to the double sum

$$\mathcal{Z}(G | \xi, \zeta) = \sum_{x \in \mathbb{N}_n} \sum_{y \in \mathbb{N}_m} \xi^x \zeta^y \mu(x, y).$$

Let $\Lambda_y = \{x \in \mathbb{N}_n \mid x \equiv 0 \pmod 2, x \le 2y\}$ be the set containing the support $\Upsilon_y = \{x \in \Lambda_y \mid \mu(x | y) > 0\}$ of the conditional distribution $\mu_y = \mu(x | y)$ of $X | Y = y$ for $y \in \mathbb{N}_m$. The set $\Lambda_y$ consists of even integers since the sum of degrees of any graph is even. The support of the joint distribution $\mu(x, y) = \mu(x | y)\mu(y)$ is thus $\Upsilon = \bigcup_{y \in \mathbb{N}_m} \Upsilon_y$, where $\mu(y) = |\mathcal{G}_y|/|\mathcal{G}| = \binom{m}{y}/2^m$ is the marginal distribution of Y and $\mathcal{G}_y = \{g \in \mathcal{G} \mid Y(g) = y\}$.

Let $\mathbf{x} : \mathcal{G} \to \{0, 1\}^n$ be the vector-valued function with components $x_w = x_w(g) = d_w(g) \pmod 2$ for $w \in \mathbf{V}$, where $d_w(g) \in \mathbb{N}_{m-1}$ is the degree of vertex $w \in \mathbf{V}$ in the subgraph $g \in \mathcal{G}$, so that $X(g) = \sum_{w \in \mathbf{V}} x_w(g)$. Let the support $\Omega_y = \{\mathbf{x}(g) \in \{0, 1\}^n \mid g \in \mathcal{G}_y\}$ for $y \in \mathbb{N}_m$ be partitioned by the sets $\Omega_{xy} = \{\mathbf{x}(g) \in \{0, 1\}^n \mid g \in \mathcal{G}_{xy}\}$ for $x \in \mathbb{N}_n$.
### 3.2. Low-rank approximation.
Let $v_y \subseteq \Upsilon_y$ for $y \in \mathbb{N}_m$ and define the low-rank approximation

$$f_{\Phi}(x, y) \propto \begin{cases} \binom{m}{y} |\Omega_{xy}| \exp\left\{\sum_{k \in v_{y0}} \phi_{yk} t_k(x)\right\} & \text{for } (x, y) \in \Upsilon \\ 0 & \text{otherwise} \end{cases}$$

to the joint distribution $\mu(x, y)$, where $v_{y0} = v_y \setminus \{0\}$. The parameter set $\Phi \in \mathcal{R}_{v_0}^{n \times m}$, where

$$\mathcal{R}_{v_0}^{n \times m} = \times_{y \in \mathbb{N}_m} \times_{k \in v_{y0}} \begin{cases} \mathbb{R} & \text{if } k \in v_{y0} \\ \{0\} & \text{otherwise} \end{cases}$$

has components $\phi_{yk} \in \mathbb{R}$ for $k \in v_{y0}$ and $\phi_{yk} = 0$ otherwise for $y \in \mathbb{N}_m$.

Consider the constrained optimization problem

$$\begin{array}{ll} \min_{\Phi \in \mathcal{R}_{v_0}^{n \times m}} & \mathrm{KL}(\mu \| f_{\Phi}) \\ \text{subject to} & \sum_{y \in \mathbb{N}_m} f_{\Phi}(0, y) = \sum_{y \in \mathbb{N}_m} \mu(0, y), \end{array}$$

where $\mu = \mu(x, y)$ and the constraint conserves the probability that a random spanning subgraph of $G$ is a cycle, i.e., is Eulerian. In particular, let $C = \{g \in \mathcal{G} | X = 0\}$ be the set of spanning subgraphs of $G$ that are cycles. The set $C \subseteq \mathcal{G}$ has cardinality $|C| = 2^{m-n+1}$ [7] and, thus,

$$\sum_{y \in \mathbb{N}_m} \mu(0, y) = \frac{1}{2^{n-1}}.$$

The constraint may be represented explicitly through the sufficient statistic $t_*(x) = 1\{x = 0\}$ in the approximating distribution

$$(4) \quad f_{\Phi_*(x, y)} \propto \begin{cases} \binom{m}{y} |\Omega_{xy}| \exp\left\{\phi_* t_*(x) + \sum_{k \in v_{y0}} \phi_{yk} t_k(x)\right\} & \text{for } (x, y) \in \Upsilon \\ 0 & \text{otherwise,} \end{cases}$$

where $\phi_* \in \mathbb{R}$ and $\Phi_* = \Phi \cup \{\phi_*\} \in \mathcal{R}_{v_0}^{n \times m} \times \mathbb{R}$. In particular, define

$$\hat{\Phi}_* = \arg \min_{\Phi_* \in \mathcal{R}_{v_0}^{n \times m} \times \mathbb{R}} \mathrm{KL}(\mu \| f_{\Phi_*})$$

as the stationary point of the convex function

$$\mathrm{KL}(\mu \| f_{\Phi_*}) = \sum_{y \in \mathbb{N}_m} \mathrm{KL}(\mu_y \| f_{\phi_{y*}}) \mu(y),$$

in which $\mu_y = \mu(x | y)$ and

$$f_{\phi_{y*}}(x | y) \propto \begin{cases} |\Omega_{xy}| \exp\left\{\phi_{*} t_{*}(x) + \sum_{k \in v_{y0}} \phi_{yk} t_{k}(x)\right\} & \text{for } x \in v_{y} \\ 0 & \text{otherwise,} \end{cases}$$

where $\phi_y = \Phi_{\cdot y}$ and $\phi_{y*} = (\phi'_y \ \phi_*)' \in \mathcal{R}_{v_{y0}}^n \times \mathbb{R}$ for $y \in \mathbb{N}_m$.

As for an exact representation of $\mu$, define

$$\hat{\Theta}_* = \arg \min_{\Theta_* \in \mathcal{R}_{\Upsilon_0}^{n \times m}} \mathrm{KL}(\mu \| f_{\Theta_*})$$
as the stationary point of the convex function

$$KL(\mu \| \hat{\Theta}_{*}) = \sum_{y \in \mathbb{N}_m} KL(\mu_y \| \hat{\theta}_{y*}) \mu(y),$$

where $\theta_y = \Theta_y$ and $\theta_{y*} = (\theta'_y \ \theta'_y)' \in \mathcal{R}_{\Upsilon_{y0}}^n \times \mathbb{R}$ for $y \in \mathbb{N}_m$, $\Upsilon_0 = \cup_{y \in \mathbb{N}_m} \Upsilon_{y0}$ and $\Upsilon_{y0} = \Upsilon_y \setminus \{0\}$ for $y \in \mathbb{N}_m$, and

$$\hat{f}_{\theta_{y*}}(x|y) \propto \begin{cases} |\Omega_{xy}| \exp\{\theta_{*}t_{*}(x) + \sum_{k \in \Upsilon_{y0}} \theta_{yk}t_{k}(x)\} & \text{for } x \in \Upsilon_{y} \\ 0 & \text{otherwise.} \end{cases}$$

The increase in information entropy that results from the $|v_0|$-rank approximation $\hat{f}_{\hat{\Theta}_{*}}$ to the probability distribution $\hat{f}_{\hat{\Theta}_{*}}$ is

$$KL(\hat{f}_{\Theta_{*}} \| \hat{f}_{\Phi_{*}}) = \sum_{y \in \mathbb{N}_m} (\phi_y - \theta_y) \cdot (E_{\hat{f}_{\phi_{y*}}} \mathbf{t} - E_{\hat{f}_{\theta_{y*}}} \mathbf{t}) \cdot \mu(y) + O\left(\sum_{y \in \mathbb{N}_m} (\|\phi_y - \theta_y\|^2 + (\phi_{*} - \theta_{*})^2) \mu(y)\right),$$

where the positive constant associated with the big-O term is, once again, guaranteed by a positive definite Hessian matrix.

### 3.3. Population moments. A solution to the equations

$$0 = \nabla_{\phi_{y*}} KL(\mu_y \| \hat{f}_{\phi_{y*}}) \Big|_{v_{y0}} = (E_{\hat{f}_{\phi_{y*}}} \mathbf{t} - E_{\mu_y} \mathbf{t}) \Big|_{v_{y0}}$$

requires evaluation of the population moments

$$E_{\mu}\left(\binom{X}{k} \mid Y=y\right) = E_{\nu}\left(\mathfrak{c}_{k}(\mathbf{x}) \mid Y=y\right) = \sum_{\substack{W \subseteq \mathbf{V} \\ |W|=k}} E_{\nu}(x_W \mid Y=y)$$

for $k \in v_{y0}$ and $y \in \mathbb{N}_m$, where $x_W = \prod_{w \in W} x_w$ and $\nu = \nu(\mathbf{x}, y)$ is the joint distribution of $\mathbf{x} \in \{0, 1\}^n$ and $Y \in \mathbb{N}_m$ over $g \in \mathcal{G}$.

Let $\mathbf{E}_W = \{(u, v) \in \mathbf{E} \mid u \in W \wedge v \in W\}$ for $W \subseteq \mathbf{V}$ and define $\mathcal{G}_W = \{g = (\mathbf{V}, e) \mid e \subseteq \mathbf{E}_W\}$. The expected value

$$E_{\nu}(x_W | Y=y) = \sum_{g_W \in \mathcal{G}_W} E_{\nu}(x_W | Y=y, g_W) \frac{\binom{|\mathbf{E}_W|}{|e_W|} \binom{m-|\mathbf{E}_W|}{y-|e_W|}}{\binom{m}{y}},$$

where $e_W \subseteq \mathbf{E}_W$ is the edge set of the spanning subgraph $g_W = (\mathbf{V}, e_W)$. Let $\mathbf{E}_{\bar{W}} = \{(u, v) \in \mathbf{E} \mid (u \in W \wedge v \notin W) \lor (u \notin W \wedge v \in W)\}$ and $\mathbf{E}_{W^c} = \{(u, v) \in \mathbf{E} \mid u \notin W \wedge v \notin W\}$ and, given the configuration of $g_W$, suppose that $|e_{\bar{W}}| \in \mathbb{N}_{y-|e_W|}$ of the $y - |e_W|$ edges remaining are allocated to $\mathbf{E}_{\bar{W}}$. Thus,

$$E_{\nu}(x_W | Y=y, g_W) = \sum_{|e_{\bar{W}}|=0}^{y-|e_W|} E_{\nu}(x_W | Y=y, g_W, |e_{\bar{W}}|) \frac{\binom{|\mathbf{E}_{\bar{W}}|}{|e_{\bar{W}}|} \binom{|\mathbf{E}_{W^c}|}{y-|e_W|-|e_{\bar{W}}|}}{\binom{m-|\mathbf{E}_W|}{y-|e_W|}}$$

since $m = |\mathbf{E}_W| + |\mathbf{E}_{\bar{W}}| + |\mathbf{E}_{W^c}|$.

Let $d_w: \mathbf{E} \mapsto \mathbb{N}_{m-1}$ be the function that returns the degree $d_w(e)$ of vertex $w \in W$ in the spanning subgraph $g = (\mathbf{V}, e) \in \mathcal{G}$ and define $\mathbb{I}_w = \{1, 3, \dots, d_w(\mathbf{E}_{\bar{W}}) - 1\}$ if $d_w(\mathbf{E}_{\bar{W}}) \equiv 0 \pmod{2}$ and $\mathbb{I}_w = \{1, 3, \dots, d_w(\mathbf{E}_{\bar{W}})\}$ if $d_w(\mathbf{E}_{\bar{W}}) \equiv 1 \pmod{2}$ for
$w \in W$ and $\mathbb{I}_W = \times_{w \in W} \mathbb{I}_w$. Let $p_w(e_W) = d_w(e_W) \pmod 2$ be the parity of vertex $w \in W$ in the spanning subgraph $g_W = (\mathbf{V}, e_W) \in \mathcal{G}_W$ so that

$$E_{\nu}(x_W | Y = y, g_W, |e_{\bar{W}}|) = \sum_{\substack{\mathbf{\iota} \in \mathbb{I}_W \\ \|\mathbf{\iota}\|_1 = |e_{W^c}|}} \prod_{w \in W} \binom{d_w(\mathbf{E}_{\bar{W}})}{\iota_w - p_w(e_W)},$$

where $|e_{W^c}| = y - |e_W| - |e_{\bar{W}}|$ and $\iota_w$ is the $w$th component of $\mathbf{\iota} \in \mathbb{I}_W$.

3.4. **Approximation error.** Let $\hat{Z}(G | \alpha, \beta)$ be the approximation

$$\hat{Z}(G | \alpha, \beta) = (2 \cosh \alpha)^n (2 \cosh \beta)^m \hat{Z}(G | \xi, \zeta)$$

to the partition function $Z$ in Equation (3) where, $\xi = \tanh \alpha$ and $\zeta = \tanh \beta$ as before, and

$$\hat{Z}(G | \xi, \zeta) = \sum_{x=0}^{n} \sum_{y=0}^{m} \xi^x \zeta^y \mathfrak{f}_{\hat{\Phi}_*}(x, y),$$

where $\mathfrak{f}_{\hat{\Phi}_*}$ is the low-rank approximation to $\mu$ in Equation (4). The relative error $R = |Z - \hat{Z}|/Z$ of an $v_0$-order approximation is thus bounded as

$$R(v_0 | \alpha, \beta) \le \left( \frac{1 - |\xi|^{n+1}}{1 - |\xi|} \right) \left( \frac{1 - |\zeta|^{m+1}}{1 - |\zeta|} \right) \frac{d_{TV}(\mu, \mathfrak{f}_{\hat{\Phi}_*})}{\mathcal{Z}(G | \xi, \zeta)},$$

where $d_{TV}$ is the total deviation metric which, itself, is bounded

$$d_{TV}(\mu, \mathfrak{f}_{\Phi_*}) \le \sqrt{\frac{1}{2} \mathrm{KL}(\mu \| \mathfrak{f}_{\Phi_*})}$$

according to Pinsker's inequality. [25]

4. DISCUSSION

In the general context of Section 2, the fundamental relation

$$|\Omega_x| \exp \left\{ \sum_{k \in \Upsilon_0} \theta_k \binom{x}{k} \right\} = \sum_{\mathbf{x} \in \Omega_x} \exp \left\{ \sum_{k \in \Upsilon_0} \theta_k \mathbf{e}_k(\mathbf{x}) \right\}$$

between $Q(\mathbf{x})$ and $P(x)$ for $\mathbf{x} \in \Omega$ and $x \in \mathbb{N}_n$ such that $X = x$ elicits the probability distribution $q_{\theta} : \{0, 1\}^n \mapsto \mathbb{R}_+ \cup \{0\}$ with expression

$$q_{\theta}(\mathbf{x}) = \begin{cases} \exp\{\sum_{k \in \Upsilon_0} \theta_k \mathfrak{e}_k(\mathbf{x}) - \Psi_q(\theta)\} & \text{for } \mathbf{x} \in \Omega \\ 0 & \text{otherwise} \end{cases}$$

that is invariant to permutations of the components of $\mathbf{x} \in \Omega$ and has a tractable cumulant function $\Psi_q(\theta) = \log \sum_{x \in \Upsilon} \exp\{\sum_{k \in \Upsilon_0} \theta_k \binom{x}{k}\}$.

As mentioned at the outset, if $Q$ is ergodic then $P'(x) \propto \binom{n}{x}$ and

$$h_{\theta}(x) \propto \binom{n}{x} \exp \left\{ \sum_{k=1}^{n} \theta_k \binom{x}{k} \right\}$$

for $x \in \mathbb{N}_n = \{0, 1, \dots, n\}$. Although there is no restriction that $Q$ itself be symmetric in the components of $\mathbf{x} \in \Omega$, the symmetric model $q_{\theta}(\mathbf{x})$ of a two-state system is sufficient for the representation of the probability distribution $P$ of the sum $X = x_1 + \dots + x_n$ of the components of $\mathbf{x}$ as a member of the family of probability distributions $h_{\theta} : \Omega \mapsto \mathbb{R}_+$ indexed by $\theta \in \mathbb{R}^{|\Upsilon_0|}$. The family $h_{\theta}$ is contained
within the exponential family of probability distributions and, since the collection $\mathcal{E} = \{\mathbf{e}_k : \{0,1\}^n \to \mathbb{N}_+ | k \in \mathbb{N}_n\}$ of elementary symmetric functions is a basis for the set $\mathcal{B} = \{\mathbf{B} : \{0,1\}^n \to \mathbb{R}\}$ of pseudo-Boolean functions, the family $h_{\theta}$ is complete.

Under the precept that probability distributions encode incomplete information about a system, we circumvented the usual obstacles encountered by classical methods of approximation by maximally sharing information from the probability distribution $Q(\mathbf{x})$ over $\Omega \subseteq \{0,1\}^n$ pertinent to the probability distribution $P(x)$ over $\Upsilon \subseteq \mathbb{N}_n$ through the statistics $\mathbf{e}_k(\mathbf{x}) = \binom{X}{k}$ for $k \in \mathbb{N}_n$.

In the case that the two-state system is ergodic, i.e., $\Omega = \{0,1\}^n$, then

The ancillary information $A_Q = \text{KL}(Q\|q_{\theta})$, where $\theta \in \mathbb{R}_{\Upsilon_0}^n$ was defined in Equation (1) as that which identified $h_{\theta} = P$, is zero if and only if $Q(\mathbf{x})$ is symmetric in the components of $\mathbf{x}$, in which case the components $x_1, \dots, x_n$ of $\mathbf{x}$ are exchangeable. As such, any statistical quantity derived by $Q$ is exclusively a function of the sum $X = x_1 + \dots + x_n$. The assumption of exchangeability would appear to be tenable in two-state physical systems with components indistinguishable except possibly for their patterns of interactivity. The probability distribution of any statistic of the configuration $\mathbf{x} \in \Omega$ is accordingly induced by the probability distribution $P$ in this context.

## REFERENCES

[1] A. D. Barbour, Lars Holst, and Svante Janson. *Poisson Approximation*. Oxford University Press, 1992.

[2] A.D. Barbour and L.H.Y. Chen. *An Introduction to Stein's Method*. Lecture Notes Series - National University of Singapore, Institute for Mathematical Sciences, Singapore University Press, 2005.

[3] A.D. Barbour and P. Hall. On the rate of poisson convergence. *Mathematical Proceedings of the Cambridge Philosophical Society*, 95:473—480, 1984.

[4] Stephen G. Brush. History of the Lenz-Ising model. *Reviews of Modern Physics*, 39:883–893, 1967.

[5] Louis H. Y. Chen. Poisson approximation for dependent trials. *Annals of Probability*, 3:534–545, June 1975.

[6] Amir Dembo and Andrea Montanari. Ising models on locally tree-like graphs. *The Annals of Applied Probability*, 20(2):565–592, 2010.

[7] Richard Diestel. *Graph Theory*. Springer Publishing Company, Incorporated, 5th edition, 2017.

[8] B. Efron. The geometry of exponential families. *The Annals of Statistics*, 6(2):362–376, 1978.

[9] Werner Ehm. Binomial approximation to the poisson binomial distribution. *Statistics abd Probability Letters*, 11(1):7–16, 1991.

[10] Michael E. Fisher. On the dimer solution of planar ising models. *Journal of Mathematical Physics*, 7(10):1776–1781, 10 1966.

[11] Andreas Galanis, Leslie A. Goldberg, and Andres Herrera-Poyatos. The complexity of approximating the complex-valued Ising model on bounded degree graphs. *SIAM Journal on Discrete Mathematics*, 36(3):2159–2204, 2022.

[12] Elad Ganmor, Ronen Segev, and Elad Schneidman. Sparse low-order interaction network underlies a highly correlated and learnable neural population
code. *Proceedings of the National Academy of Sciences*, 108(23):9679-9684, 2011.

[13] I. J. Good. Maximum entropy for hypothesis formulation, especially for multidimensional contingency tables. *The Annals of Mathematical Statistics*, 34(3):911-934, 1963.

[14] J. J. Hopfield. Neural networks and physical systems with emergent collective computational abilities. *Proceedings of the National Academy of Sciences of the United States of America*, 79:2554-2558, 1982.

[15] Ernst Ising. Beitrag zur theorie des ferromagnetismus. *Zeitschrift für Physik*, 31:253-258, 1925.

[16] E. T. Jaynes. Information theory and statistical mechanics. *Physical Review*, 106(4):620-630, 1957.

[17] E. T. Jaynes. Information theory and statistical mechanics. ii. *Physical Review*, 108(2):171-190, 1957.

[18] Edwin T. Jaynes. *Probability Theory: The Logic of Science*. Cambridge university press, 2003.

[19] Mark Jerrum and Alistair Sinclair. Polynomial-time approximation algorithms for the ising model. *SIAM Journal on Computing*, 22(5):1087-1116, 1993.

[20] M. Kac and J. C. Ward. A combinatorial solution of the two-dimensional ising model. *Phys. Rev.*, 88:1332-1337, Dec 1952.

[21] P. W. Kasteleyn. Dimer statistics and phase transitions. *Journal of Mathematical Physics*, 4(2):287-293, 02 1963.

[22] Bruria Kaufman. Crystal statistics. ii. partition function evaluated by spinor analysis. *Phys. Rev.*, 76:1232-1243, Oct 1949.

[23] Andrei Nikolaevich Kolmogorov. *Foundations of the Theory of Probability*. 1933.

[24] S. Kullback. *Information Theory and Statistics*. Wiley publication in mathematical statistics. Wiley, 1959.

[25] S. Kullback. A lower bound for discrimination information in terms of variation. *IEEE Trans. on Information Theory*, 13(1):126-127, 1967.

[26] S. Kullback and R. A. Leibler. On information and sufficiency. *The Annals of Mathematical Statistics*, 22(1):79-86, 1951.

[27] Lucien Le Cam. An approximation theorem for the Poisson binomial distribution. 10:1181-1197, 1960.

[28] Wilhelm Lenz. Beiträge zum verständnis der magnetischen eigenschaften in festen körpern. *Physikalische Zeitschrift*, 21:613-615, 1920.

[29] Jingcheng Liu, Alistair Sinclair, and Piyush Srivastava. Fisher zeros and correlation decay in the ising model. *Journal of Mathematical Physics*, 60(10):103304, 10 2019.

[30] Michael W. Macy, Boleslaw K.j Szymanski, and Janusz A. Hołyst. The ising model celebrates a century of interdisciplinary contributions. *npj Complexity*, 1:10, 2024.

[31] Gordon F. Newell and Elliott W. Montroll. On the theory of the ising model of ferromagnetism. *Reviews of Modern Physics*, 25:353-389, 1953.

[32] Charles M. Newman. Zeros of the partition function for generalized ising systems. *Communications on Pure and Applied Mathematics*, 27(2):143-159, 1974.
[33] Lars Onsager. Crystal statistics. i. a two-dimensional model with an order-disorder transition. *Physical Review*, 65:117-149, 1944.

[34] H. Peierls. Ising's model of ferromagnetism. *Proceedings of Cambridge Philosophical Society*, 32:477-481, 1936.

[35] Erol A. Peköz, Adrian Röllin, Vydas Čekanavičius, and Michael Shwartz. A three-parameter binomial approximation. *Journal of Applied Probability*, 46(4):1073-1085, 2009.

[36] Han Peters and Guus Regts. Location of zeros for the partition function of the Ising model on bounded degree graphs. *Journal of the London Mathematical Society*, 101(2):765-785, 2020.

[37] Elad Schneidman, Michael J. Berry, Ronen Segev, and William Bialek. Weak pairwise correlations imply strongly correlated network states in a neural population. *Nature*, 440(7087):1007-1012, 2006.

[38] H. S. Seung, H. Sompolinsky, and N. Tishby. Statistical mechanics of learning from examples. *Physical Review A*, 45:6056-6091, Apr 1992.

[39] C. E. Shannon. A mathematical theory of communication. *The Bell System Technical Journal*, 27(3):379-423, 1948.

[40] Alistair Sinclair, Piyush Srivastava, and Marc Thurley. Approximation algorithms for two-state anti-ferromagnetic spin systems on bounded degree graphs. In *Proceedings of the Twenty-Third Annual ACM-SIAM Symposium on Discrete Algorithms*, SODA '12, page 941-953, USA, 2012. Society for Industrial and Applied Mathematics.

[41] Alistair Sinclair, Piyush Srivastava, and Marc Thurley. Approximation algorithms for two-state anti-ferromagnetic spin systems on bounded degree graphs. *Journal of Statistical Physics*, 155:666-686, 2014.

[42] Allan Sly and Nike Sun. Counting in two-spin models on d-regular graphs. *The Annals of Probability*, 42(6):2383-2416, 2014.

[43] B. L. van der Waerden. Die lange Reichweite der regelmäßigen Atomanordnung in Mischkristallen. *Zeitschrift fur Physik*, 118(7-8):473-488, July 1941.

[44] Martin J. Wainwright and Michael I. Jordan. Graphical models, exponential families, and variational inference. *Found. Trends Mach. Learn.*, 1(1-2):1-305, January 2008.

[45] D. J. A. Welsh. *The computational complexity of some classical problems from statistical physics*, pages 307-321. Oxford University Press, 1990.

[46] Cristian Zanoci, Nima Dehghani, and Max Tegmark. Ensemble inhibition and excitation in the human cortex: An ising-model analysis with uncertainties. *Phys. Rev. E*, 99:032408, Mar 2019.