\documentclass{jors}

%% Set the header information
\pagestyle{fancy}
\definecolor{mygray}{gray}{0.6}
\renewcommand\headrule{}
\rhead{\footnotesize 3}
\rhead{\textcolor{gray}{UP JORS software Latex paper template version 0.1}}

\usepackage{graphics,color,graphicx,amsmath,hyperref}
\usepackage{xr,natbib}
\bibliographystyle{vancouver}
\usepackage{amsbsy}

\newcommand{\norm}[1]{\left\Vert#1\right\Vert}
\newcommand{\abs}[1]{\left\vert#1\right\vert}
\newcommand{\set}[1]{\left\{#1\right\}}
\newcommand{\eps}{\varepsilon}
%\bibliographystyle{abbrvnat}
\graphicspath{{images}}

\begin{document}

% COMMANDS -------------------------------------------------------
\newcommand{\mr}[1]{\mathrm{#1}}
\newcommand{\mb}[1]{\mathbf{#1}}
\newcommand{\br}[1]{\left<#1\right>}
\newcommand{\bl}[1]{\left|#1\right|}
\newcommand{\mc}[1]{\mathcal{#1}}
\newcommand{\tb}[1]{\textcolor{blue}{#1}}
\newcommand{\tr}[1]{\textcolor{red}{#1}}
\newcommand{\tg}[1]{\textcolor{green}{#1}}
\newcommand{\si}[0]{{\rm s}_{\rm i}}
\newcommand{\sj}[0]{{\rm s}_{\rm j}}
\newcommand{\bs}[1]{\boldsymbol{#1}}
\newcommand{\rs}[0]{{\rm s}}
\newcommand{\rk}[0]{{\rm k}}

\rule{\textwidth}{1pt}

\section*{(1) Overview}
\vspace{0.5cm}

\section*{Convenient interface to inverse Ising (ConIII): A Python 3 package for solving Ising-type maximum entropy models}

\section*{Paper Authors}
1. Lee, Edward D.; \\
2. Daniels, Bryan C.

\section*{Paper Author Roles and Affiliations}
1. Department of Physics, 142 Sciences Dr, Cornell University, Ithaca NY 14853 \\
2. ASU--SFI Center for Biosocial Complex Systems, Arizona State University, Tempe, AZ 85287

\section*{Abstract}

ConIII (pronounced CON-ee) is an open-source Python project providing a simple interface to solving the pairwise and higher order Ising model and a base for extension to other maximum entropy models. We describe the maximum entropy problem and give an overview of the algorithms that are implemented as part of ConIII (\url{https://github.com/eltrompetero/coniii}) including Monte Carlo histogram, pseudolikelihood,  minimum probability flow, a regularized mean field method, and a cluster expansion method. Our goal is to make a variety of maximum entropy techniques accessible to those unfamiliar with the techniques and accelerate workflow for users.

\section*{Keywords}
inverse Ising, maxent, maximum entropy, inference

\section*{Introduction}
Many biological and social systems are characterized by collective behavior: the correlated pattern of neural firing \cite{Schneidman:2006he}, protein diversity in the immune system \cite{Mora:2010jxb}, conflict participation in monkeys \cite{Daniels:2017cq}, flocking in birds \cite{Bialek:2012cs}, statistics of letters in words \cite{Stephens:2010hi}, or consensus voting in the US Supreme Court \cite{Lee:2015ev,Lee:ky}. Statistical physics is a natural approach to probing such systems precisely because they are collective \cite{Castellano:2009ce}.
Recently, the development of numerical, analytic, and computational tools have made it feasible in these large collective systems to solve for the maximum entropy (maxent) model that reproduces system behavior, corresponding to solving an ``inverse problem.''
This approach contrasts with the typical problem in statistical physics where one postulates the microscopic model (the Hamiltonian) and works out the physical behavior of the system. In the inverse problem, we  find the parameters that correspond to observed behavior of a known system. In many cases, this is a very difficult problem to solve and does not have an analytical solution, and we must rely on analytic approximation and numerical techniques to estimate the parameters.

The pairwise maxent model, the Ising model, has been of particular interest because of its simplicity and generality. A variety of algorithms have been proposed to solve the inverse Ising problem, but different approaches are disparately available on separate code bases in different coding languages, which makes comparison difficult and pedagogy more complicated. With ConIII, it is possible to solve the inverse Ising problem with a variety of algorithms in just a few lines of code. 

ConIII is intended to provide a centralized resource for the inverse Ising problem and easy extension to other maxent problems. Although some of the implemented algorithms are specific to the pairwise Ising model, maxent models with arbitrary combinations of higher order constraints can be solved as well by specifying the particular constraints of the maxent model of interest. 
%The modular structure of ConIII was inspired by the machine learning package scikit-learn, where each algorithm is implemented in a class.
%There are a variety of techniques for approaching the inverse problem|often in context of a particular model|but they are often presented as different approaches without clear relations to one another. Instead of focusing on the details of each technique separately, we relate these techniques to one another and point out when one would fail or be difficult to apply.

%We hope to provide a (soft) introduction to the analytic, numerical, and computational techniques used to solve these problems to make them accessible to students especially at the graduate level in physics and beyond (and even advanced undergraduates). With the goal of furthering general understanding, we aim not to be rigorous but to be accurate and comprehensible, pointing to the relevant literature whenever possible. Although we strive to be as clear and explicit as possible across this guide, a good background in mathematics or physics will be extremely helpful. In some places, we have left out steps in the derivations, and we encourage the reader to work out the missing steps.

In the first few sections of this paper, we give a brief overview of maxent and describe at a high level the algorithms implemented in this package. For those unfamiliar with maxent, we also provide some useful references like Ref~\cite{Bialek:2012ueb} and the appendix of \cite{Lee:2015ev}. For those seeking more detail about the implemented algorithms, we provide references specific to each algorithm section. Then, we describe the architecture of the package and how to contribute.

\section*{What is maximum entropy?}
Shannon introduced the concept of information entropy in his seminal paper about communication over a noisy channel \cite{Shannon:1948wk}. Information entropy is the unique measure of uncertainty that follows from insisting on elementary principles of consistency. According to Shannon, the entropy over the probability distribution $p(\rs)$ of possible discrete configurations $\mathcal S$ of a system is
\begin{align}
	S[p] &= -\sum_{{\rm s}\in \mathcal{S}} p({\rm s}) \log p({\rm s}).
\end{align}
These configurations could be on-off patterns of firing in neurons, the arrangement of letters in a word, or the orientation of spins in a material.

When there is no structure in the distribution, meaning that the probability is uniform, entropy is at a maximum. In the context of communication theory as Shannon first discussed, this means that there is no structure to exploit to make a prediction about the next part of an incoming message; thus, maximum entropy means that each new part of the message is maximally ``surprising.'' At the other extreme, when the message consists of the same bit over and over again, we can always guess at the following part of the message and the signal has zero entropy. In the context of modeling, we use entropy not to refer to the difficulty of the message, but to our state of knowledge about it. Entropy precisely measures our uncertainty about the configuration in which we expect to find the system.

Maximum entropy, or maxent, is the formal framework for building models that are consistent with statistics from the data but otherwise as structureless as possible \cite{Bretthorst:2003ua,Jaynes:1957fy}.
We begin by choosing a set of $K$ useful or important features from the data $f_{\rm k}(\rs)$ that should be true for the model that we are trying to build. These could be whether or not a set of neurons fire together in a temporal bin or the pairwise coincidence for primates in a conflicts. The average of this feature across the data set $\mathcal{D}$ with $R$ samples is
\begin{align}
	\br{f_{\rm k}}_{\rm data} &= \frac{1}{R}\sum_{{\rm s}\in \mathcal{D}} f_{\rm k}({\rm s}).
\end{align}
According to the model in which each observation s occurs with some probability $p(\rm s)$, the same average is calculated over all possible states
\begin{align}
	\br{f_{\rm k}} &= \sum_{\rm s\in\mathcal{S}} p(\rs)f_{\rm k}(\rs).
\end{align}

We assert that the model should fit the $K$ features while maximizing entropy. The standard procedure is to solve this by the method of Langrangian multipliers. We construct the Langrangian functional $\mathcal{L}$ by introducing the multipliers $\lambda_{\rm k}$.
\begin{align}
	\mathcal{L}[p] &= -\sum_{\rm s\in\mathcal{S}} p(\rs)\log p(\rs) - \sum_{\rm k=1}^{K} \lambda_{\rm k} \big(\br{f_{\rm k}}-\br{f_{\rm k}}_{\rm data}\big)\label{eq:Lagrang}
\end{align}


%In the notation of statistical physics, the Langrangian is the Helmholtz free energy describing the competition between entropy and the structure described in the Hamiltonian $E$ in equilibrium (See Appendix) \footnote{If we fix the average energy of the system, we get the microcanonical ensemble (See Appendix). which is equivalent to the axiom that all microstates have equal probability (that is the maxent distribution).}.
%\begin{align}
%	F &= S - \br{E}\\
%	E &= -\sum_{\rm k}^K\lambda_{\rm k}f_{\rm k}(s)
%\end{align}
%This formulation makes clear the fundamental connection that statistical mechanics is an inference procedure using the maximum entropy principle \cite{Jaynes:1957fy}.

Then, we solve for the fixed point by taking the derivative with respect to $\lambda_{\rm k}$. 
The resulting maxent model is a Boltzmann distribution over states:
\begin{align}
    \label{eq:energy}
	p(\rs) &= \left.e^{-E(\rs)}\right/Z,\\
\intertext{with relative negative log-likelihood (also known as the energy or Hamiltonian)}
E(\rs) &= -\sum_{\rm k=1}^K\lambda_{\rm k}f_{\rm k}(\rs),
\intertext{and normalization factor (also known as the partition function)}
	Z &= \sum_{\rm s\in\mathcal{S}} e^{-E(\rs)}.
\end{align}
Entropy is a convex function of $p$ and the constraints are linear with respect to $p$, so the problem is convex and the maxent distribution unique.
Readers familiar with statistical physics will recognize this as an alternative derivation of the microcanonical ensemble, demonstrating that statistical mechanics can be viewed as an inference procedure using the maxent principle \cite{Jaynes:1957fy}.

Finding the parameters $\lambda_{\rm k}$ that match the constraints $\br{f_{\rm k}}_{\rm data}$ is equivalent to minimizing the Kullback-Leibler divergence between the model and the data \cite{Cover:2006tl}
\begin{align}
	D_{\rm KL}(p_{\rm data}||p) &= \sum_\rs p_{\rm data} \log\left(\frac{p_{\rm data}(\rs)}{p(\rs)}\right)\\
	\frac{\partial D_{\rm KL}}{\partial \lambda_{\rm k}} &= \sum_{\rs} p_{\rm data}(\rs) \frac{\partial (-E(\rs)-\log Z)}{\partial \lambda_{\rm k}} = 0 \nonumber \\
	\implies  \br{f_{\rm k}}_{\rm data} &= \br{f_{\rm k}}\label{eq:fk-fk}.
\end{align}
 In other words, the parameters of the maxent model are the ones that minimize the information theoretic ``distance'' to the distribution of the data given the constraints. Note that these parameters are given by the data: once the constraints have been chosen, there is a single maxent solution, with no free parameters.

\subsection*{The Ising model}
The Ising model is a statistical physics model of magnetism \cite{Ising:1924vf}. It consists of a set of spins $\{\rs_{\rm i}\}$ with 2 possible orientations (up and down), each responds to its own external magnetic field $h_{\rm i}$ and each pair is coupled to each other with pairwise coupling $J_{\rm ij}$. The strength of the magnetic field determines the tendency of each of the spins to orient in a particular direction and the couplings determine whether the spins tend to point together ($J_{\rm ij}>0$) or against each other ($J_{\rm ij}<0$). Typically, neighbors are defined as spins that interact with one another given by some underlying network structure. Figure \ref{gr:ising} shows a fully-connected example.

\begin{figure}[tb]\centering
	\includegraphics[width=.65\linewidth,clip,trim={100 70 70 60}]{ising_example.pdf}
\caption{Example of a fully connected pairwise Ising model with positive and negative couplings. Each spin $\rs_{\rm i}$ (circle) can take one of two states (black or white, corresponding to $-1$ and $1$) and is connected to every other spin in the system with a positive (red) or negative (blue) coupling. These states could describe the on-off patterns of firing in neurons, the orientation of spins in a material, or if each spin is no longer binary the arrangement of letters in a word (a Potts model).}
\label{gr:ising}
\end{figure}

The energy of each configuration determines its probability via Eq \eqref{eq:energy},
\begin{align}
	E(\rs) &= -\sum_{\rm i<j}^N J_{\rm ij}\si\sj -\sum_{\rm i=1}^Nh_{\rm i}\si,
\end{align}
such that lower energy states are more probable.
%States with lower energy are more likely so spin tend to orient along the direction of the local field $h_{\rm i}^{\rm local} = \sum_{\rm j}J_{\rm ij}\sj +h_{\rm i}$. 

We can derive the Ising model from the perspective of maxent. Fixing the the means and pairwise correlations to those observed in the data
\begin{align}
	\br{\si} &= \br{\si}_{\rm data} \\
	\br{\si\sj} &= \br{\si\sj}_{\rm data}
\end{align}
we go through the procedure of constructing the Langrangian from Eq~\ref{eq:Lagrang}
\begin{align}
	\mathcal{L}[p] &= -\sum_\rs p(\rs)\log p(\rs) +\sum_{\rm i<j}^N J_{\rm ij}\br{\si\sj} +\sum_{\rm i}^Nh_{\rm i}\br{\si}\\
	\frac{\partial\mathcal{L}[p]}{\partial p(\rs)} &= -\log p(\rs)-1 +\sum_{\rm i<j}^N J_{\rm ij}\si\sj +\sum_{\rm i}^Nh_{\rm i}\si\\
	\log p(\rs) &= -1 +\sum_{\rm i<j}^N J_{\rm ij}\si\sj +\sum_{\rm i}^Nh_{\rm i}\si\label{eq:E-1}\\
	p(\rs) &= \left.e^{-E(\rs)}\right/Z\\
\intertext{where the $-1$ in Eq \ref{eq:E-1} has been absorbed into the normalization factor}
	Z &= \sum_{\rs} e^{-E(\rs)}.
\end{align}
such that the probability distribution is normalized $\sum_{\rs} p(s)=1$.
Thus, the resulting model is exactly the Ising model mentioned earlier.

Despite the simplicity of the Ising model, the structure imposed by the discrete nature of the spins means that finding the parameters is challenging analytically and computationally. In the last few years, numerous techniques have been suggested for solving the inverse Ising problem exactly or approximately \cite{Nguyen:2017ww}. We have implemented some of them in ConIII and designed a package structure to make it easily extensible to include more methods.  Here, we briefly describe the algorithms that are part of the first official version of the package. The goal is to give the user a sense of how they work without getting bogged down in heavy detail. For more detail, we suggest perusing the papers referenced in each section or the review \cite{Nguyen:2017ww}. For a complete beginner, it may be useful to first get familiar with a slower introduction like in the Appendices of Ref \cite{Lee:2015ev}, Ref \cite{Bialek:2012ueb}, or Ref \cite{Bretthorst:2003ua}.



\section*{Inverse Ising methods implemented in ConIII}
\subsection*{Enumeration}
The na\"{i}ve approach that only works for small systems is to write out the equations from Eq~\ref{eq:fk-fk} and solve them numerically. After writing out all $K$ equations,
\begin{equation}
	\br{f_{\rm k}} = -\frac{\partial \ln Z}{\partial \lambda_{\rm k}} = \br{f_{\rm k}}_{\rm data} \\,
\end{equation}
we can use any standard root-finding algorithm to find the parameters $\lambda_{\rm k}$.
This approach, however, involves enumerating all states of the system, whose number grows exponentially with system size.
%[Simple example with logical gates?]

For the Ising model, writing down the equations has a number of steps $\mathcal{O}(K^2 2^{N})$, where $K$ is the number of constraints and $N$ the number of spins.  Each evaluation of the objective in the root-finding algorithm will be of the same order. For relatively small systems, around $N\leq15$, this approach is feasible on a typical desktop computer and is a good way to test the results of a more complicated algorithm.

This approach is implemented by the {\tt Enumerate} class.


\subsection*{Monte Carlo Histogram (MCH)}
Perhaps the most straightforward and most expensive computational approach is Monte Carlo Markov Chain (MCMC) sampling.  A series of states sampled from a proposed $p(\rs)$ is produced by MCMC to approximate $\br{f_{\rm k}}$ and determine how close we are to matching $\br{f_{\rm k}}_{\rm data}$.  The parameters are then adjusted using a learning rule, and both sampling and learning are repeated until a stopping criterion is met.
%The magic behind MCMC approaches is that only relative probabilities of states need to be known to return a sample of the distribution. Thus, we do not need to calculate the partition function, but can instead compute the difference in energy between two states. 
This can be combined with a variety of approximate gradient descent methods to reduce the number of sampling steps by predicting how the distribution will change if we modify the parameters slightly.
The particular technique implemented in ConIII is the Monte Carlo Histogram (MCH) method \cite{Broderick:2007wq}.

Since the sampling step is expensive, the idea behind MCH is to reuse a sample for more than one gradient descent step \cite{Broderick:2007wq}. Given that we have a sample with probability distribution $p(\rs)$ generated with parameters $\lambda_{\rm k}$, we would like to estimate the proposed distribution $p'(\rs)$ from adjusting our parameters $\lambda_{\rm k}' = \lambda_{\rm k}+\Delta\lambda_{\rm k}$. We can leverage our current sample to make this extrapolation.
\begin{align}
	p' &= \frac{p'}{p}p\\
	p'(\rs)	&= \frac{Z}{Z'}e^{\sum_\rk \Delta\lambda_k f_\rk(\rs)} p(\rs)\\
\intertext{To estimate the average,}
	\sum_\rs p'(\rs) f_{\rm k}(\rs) &= \frac{Z}{Z'} \sum_\rs p(\rs) e^{\sum_{\rm k} \Delta\lambda_{\rm k} f_{\rm k}(\rs)} f_{\rm k}(\rs)\\
\intertext{To be explicit about the fact that we only have a sampled approximation to $p$, we replace $p$ with the sample distribution.}
	\br{f_{\rm k}}' &= \frac{Z}{Z'} \br{e^{\sum_{\rm k} \Delta\lambda_{\rm k} f_{\rm k}(\rs)} f_{\rm k}(\rs)}_{\rm sample}
\end{align}
Likewise, the ratio of the partition function can be estimated
\begin{align}
	\frac{Z}{Z'} \approx 1\left/\br{e^{\sum_{\rm k} \Delta\lambda_{\rm k} f_{\rm k}(\rs)}}_{\rm sample}\right.
\end{align}

At each step, we update the Lagrangian multipliers $\{\lambda_\rk\}$ while being careful to stay within the bounds of a reasonable extrapolation. One suggestion is to update the parameters with some inertia \cite{Tkacik:2006vq}
\begin{align}
	\Delta\lambda_\rk(t+1) &= \Delta \lambda_\rk(t) + \epsilon \Delta\lambda_\rk(t-1)\label{eq:mch learn1}\\
	\Delta \lambda_\rk(t) &= \eta\left(\br{f_\rk}'-\br{f_\rk}\right)\label{eq:mch learn2}
\end{align}
This has a fixed point at the correct parameters.
%Another heuristic is to begin with a small sample size and increase the size of the sample set as we refine our estimates of the parameters.

In practice, MCH can be difficult to tune properly and one must check in on the progress of the algorithm often. One issue is choosing how to set the learning rule parameters $\eta$ and $\epsilon$. One suggestion for $\eta$ is to shrink it as the inverse of the number of iterations \cite{Tkacik:2006vq}. Another issue is that parameters cannot be changed by too much when using the MCH approximation step or the extrapolation to $\lambda_{\rm k}'$ will be inaccurate and the algorithm will fail to converge. In ConIII, this can be controlled by setting a bound on the maximum possible change in each parameter  $\Delta\lambda_{\rm max}$ and restricting the norm of the vector of change in parameters $\sum_{\rm k} \sqrt{\Delta\lambda_{\rm k}^2}$. Another issue is setting the parameters of the MCMC sampling routine. Both the burn time (the number of iterations before starting to sample) and sampling iterations (number of iterations between samples) must be large enough that we are sampling from the equilibrium distribution.  Typically, these are found by measuring how long the energy or individual parameter values remain correlated as MCMC progresses. The parameters may need to be updated during the course of MCH because the sampling parameters may need to change with the estimated parameters of the model. For some regimes of parameter space, samples are correlated over long times and alternative sampling methods like Wolff or Swendsen-Wang would vastly reduce time to reach the equilibrium distribution although these are not included in the current release of ConIII. We do not discuss these sampling details here, but see Refs \cite{MacKay:2005wc,Newman:1999wu} for examples.

The main computational cost for MCH lies in the sampling step. For each iteration of MCH, the runtime is proportional to the number of samples $n$, number of MCMC iterations $T$, and the number of constraints for the Ising model $N^2$,
$\mathcal{O}(T n N^2)$, whereas the MCH estimate is relatively quick $\mathcal{O}(t n N^2)$ because the number of MCH approximation steps needed to converge is much smaller than the number of MCMC sampling iterations $t<<T$.
%What is the runtime for the learning rules Eqs \ref{eq:mch learn1} and \ref{eq:mch learn2}? This a much harder question that gets into optimization theory.

MCH is implemented in the {\tt MCH} class.

\subsection*{Pseudolikelihood}
The pseudolikelihood approach is an analytic approximation to the likelihood that drastically reduces the computational complexity of the problem and is exact as $N \rightarrow \infty$ \cite{Aurell:2012hi}. We calculate the conditional probability of each spin $\rs_{\rm i}$ given the rest of the system $\{\rs_{\rm j\neq i}\}$
\begin{align}
	p\left(\rs_{\rm i}|\{\rs_{\rm j\neq i}\}\right) &= \left( 1+e^{-2\rs_{\rm i} \left(h_{\rm i}+\sum_{j\neq i}J_{\rm ij}\rs_{\rm j}\right)} \right)^{-1}
\end{align}
Taking the logarithm, we define the approximate log-likelihood by summing over data points indexed by r:
\begin{align}
	f(h_{\rm i},\{J_{\rm ij}\}) &= \sum_{\rm r=1}^R \ln p\left(\left.\rs_{\rm i}^{(\rm r)}\right|\{\rs_{\rm j\neq i}\}^{(\rm r)}\right).
\end{align}
In the limit where the ensemble is well sampled, the average over the data can be replaced by an average over the ensemble:
\begin{align}
	f(h_{\rm i},\{J_{\rm ij}\}) &= \sum_{\rm s} \ln p\left(\left.\rs_{\rm i}\right|\{\rs_{\rm j\neq i}\}\right)p\left(\rs;h_{\rm i},\{J_{\rm ij}\}\right).
\end{align}

To find the point of maximum likelihood for a single spin $\rs_{\rm i}$, we calculate the analytical gradient and Hessian, $\partial f/\partial J_{\rm ij}$ and $\partial^2 f/\partial J_{\rm ij}\partial J_{\rm i'j'}$ for a Newton conjugate-gradient descent method. After maximizing likelihood for all spins, the maximum likelihood parameters may not satisfy the symmetry $J_{\rm ij}=J_{\rm ji}$. We impose the symmetry by insisting that
\begin{align}
	J_{\rm ij}' &= (J_{\rm ij}+J_{\rm ji})/2.
\end{align}

Pseudolikelihood is extremely fast and often surprisingly accurate. Each calculation of the gradient is order $\mathcal{O}(RN^2)$ and Hessian $\mathcal{O}(RN^3)$, which must be done for all $N$. With analytic forms for the gradient and Hessian, the conjugate-gradient descent method tends to converge quickly.

Pseudolikelihood for the Ising model is implemented in {\tt Pseudo}.

\subsection*{Minimum Probability Flow (MPF)}
Minimum probability flow involves analytically approximating how the probability distribution \textit{changes} as we modify the \textit{configurations} \cite{Sohl-Dickstein:2009tt,SohlDickstein:2011im}. In the methods so far mentioned, the approach has been to maximize the objective (the likelihood function) by immediately taking the derivative with respect to the parameters. With MPF, we first posit a set of dynamics that will lead the data distribution to equilibrate to that of the model. When these distributions are equivalent, then there is no ``probability flow'' between them. This technique is closely related to score matching, where we instead have a continuous state space and can directly take the derivative with respect to the states without specifying dynamics \cite{Hyvarinen:2007ed}.

First note that Monte Carlo dynamics (satisfying ergodicity and detailed balance) would lead to equilibration to the stationary distribution. One such transition matrix suggested in Ref \cite{SohlDickstein:2011im} is
\begin{align}
	\dot{p}_{\rm s} &= \sum_{\rm s'\neq s} \Gamma_{\rm ss'} p_{\rm s'} -\sum_{\rm s'\neq s} \Gamma_{\rm s's} p_{\rm s}\\
	\Gamma_{\rm ss'} &= g_{\rm ss'}\exp\left[ \frac{1}{2}\left( E_{\rm s'}-E_{\rm s} \right) \right]
\end{align}
with transition probabilities $\Gamma_{\rm ss'}$ from state ${\rm s'}$ to state s. The connectivity matrix $g_{\rm ss'}$ specifies whether there is edge between states s and ${\rm s'}$ such that probability can flow between them. By choosing a sparse $g_{\rm ss'}$ while not breaking ergodicity, we can drastically reduce the computational cost of computing this matrix.

Imagine that we start with the distribution over the states as given by the data and run the Monte Carlo dynamics. When data and model distributions are different, probability will flow between them and indicate that the parameters must be changed.  By minimizing a derivative of the Kullback-Leibler divergence, we measure how the difference between the model and the states in the data $\mathcal{D}$ changes when the dynamics are run for an infinitesimal amount of time.
\begin{align}
	L(\{\lambda_{\rm k}\}) \equiv \partial_t D_{\rm KL}(p^{(0)}||p^{(t)}\left(\{\lambda_{\rm k}\}\right)) &= \sum_{\rm s \not\in \mathcal{D}} \dot{p}_{\rm s}(\lambda_{\rm k})
\end{align}
The idea is that this derivative is also minimized with optimal parameters: the MPF algorithm looks for a minimum of the objective function $L$.

%MPF satisfies a number of useful properties:

For the Ising model, each evaluation of the objective function where $\Gamma_{\rm s s'}$ connects each data state with $G$ neighbors has runtime $\mathcal{O}(RGN^2)$. In a large fully connected system, $G\sim 2^N$ would be prohibitively large so a sparse choice is necessary.

MPF is implemented in the {\tt MPF} class.

\subsection*{Regularized mean-field method}
One attractively simple and efficient approach uses a regularized version of
mean-field theory.  In the inverse Ising problem, mean-field theory is equivalent
to treating each binary individual as instead having a continuously varying state
(corresponding to its mean value).  The inverse problem then turns into simply inverting
the correlation matrix $C$ \cite{Monasson:2011fo}:
\begin{equation}
\label{meanFieldSolution}
J^{\mathrm{mean-field}}_{\rm ij} =
    - \frac{ (C^{-1})_{\rm ij} }{ \sqrt{p_{\rm i}(1-p_{\rm i})p_{\rm j}(1-p_{\rm j})} },
\end{equation}
where
\begin{equation}
C_{\rm ij} = \frac{ p_{\rm ij} - p_{\rm i} p_{\rm j} }{ \sqrt{p_{\rm i}(1-p_{\rm i})p_{\rm j}(1-p_{\rm j})} },
\end{equation}
and where $p_{\rm i}$ corresponds to the frequency of individual i being
in the active ($+1$) state and $p_{\rm ij}$ is the frequency of the pair
i and j being simultanously in the active state.

A simple regularization scheme in this case is to discourage large values in the interaction
matrix $J_{\rm ij}$.  This corresponds to putting more weight on solutions that are closer to
the case with no interactions (independent individuals).  A particularly convenient form
adds the following term, quadratic in $J_{\rm ij}$, to the negative log-likelihood:
\begin{equation}
\gamma \sum_{\rm i} \sum_{\rm i<j} J_{\rm ij}^2 p_{\rm i} (1-p_{\rm i}) p_{\rm j} (1-p_{\rm j}).
\end{equation}
In this case, the regularized version of the mean-field solution in \eqref{meanFieldSolution}
can be solved analytically, with the slowest computational step coming from the inversion
of the correlation matrix.  For details, see Refs.~\cite{Daniels:2017cq,Barton:2013fja}.

The idea is then to vary the regularization strength $\gamma$ to move between the
non-interacting case ($\gamma \rightarrow \infty$) and the naively calculated
mean-field solution \eqref{meanFieldSolution} ($\gamma \rightarrow 0$).
While there is no guarantee that varying this one parameter will produce solutions that are
good enough to ``fit within error bars,'' this approach has been successful in at least
one case of fitting social interactions \cite{Daniels:2017cq}.

The inversion of the correlation matrix is relatively fast, bounded by $\mathcal{O}(N^3)$. Finding the optimal $\gamma$, involves Monte Carlo sampling from the model distribution, which has computational cost similar to MCH. It is, however, much more efficient because we are only optimizing a single parameter.

This is implemented in {\tt RegularizedMeanField}.



\subsection*{Cluster expansion}

Adaptive cluster expansion \cite{Monasson:2011fo,Monasson:2011fo,Barton:2013fja}
iteratively calculates terms in the
cluster expansion of the entropy $S$:
\begin{equation}
S - S_0 = \sum_\Gamma \Delta S_\Gamma,
\end{equation}
where the sum is over clusters $\Gamma$ and in the exact case
includes all $2^N - 1$ possible nonempty subsets of individuals in the system.  In the simplest version of the expansion,
one expands around $S_0 = 0$.  In some cases it can be more advantageous to
expand around
the independent individual solution or one of the mean-field solutions
described in the previous section \cite{Barton:2013fja}.

The inverse Ising problem is solved independently
on each of the clusters, which can be done exactly when the
clusters are small.  These results are used to construct a full
interaction matrix $J_{\rm ij}$.
The expansion starts with small clusters and expands to use larger
clusters, neglecting any clusters whose
contribution $\Delta S_\Gamma$ to the entropy falls below a threshold.
To find the best solution that does not overfit,
the threshold is initially set at a large value and then lowered,
gradually including more clusters in the expansion, until samples from
the resulting $J_{\rm ij}$ fit the desired statistics of the data sufficiently well.

The runtime will depend on the size of clusters included in the expansion.
If the expansion is truncated at clusters of size $n$, the worst-case runtime
would be $\mathcal{O}\left(\binom{N}{n} 2^n\right)$.  The point is that $S$ can often
be accurately estimated even when $n \ll N$.

The adaptive cluster expansion method is implemented in the {\tt ClusterExpansion} class.

%\subsection*{Bethe/Kikuchi free energy/cavity methods}
%Another approach involves a cluster-expansion of the free energy, also known as the cavity method. This has not yet been implemented in ConIII.
%The cavity method involves use of the marginalized distribution.

%\section*{Samplers}
%For estimating statistics of a given probability distribution $p(\rs)$, we have implemented in ConIII two versions of the well-known Metropolis algorithm. One is specific to the Ising model ({\tt FastMCIsing}) and the other ({\tt Metropolis}) can sample a system given an arbitrary function that calculates the energy of each state.
%
%In some parameter regimes where spins are tightly correlated, the Metropolis algorithm is very inefficient because the time to flip large numbers of spins simultaneously is extremely long. Cluster sampling like Wolff or Swendsen-Wang are much more efficient and will be included in future releases of ConIII.

%\section*{The finite sample problem}
%A fundamental problem in model inference is that uncertainty coming
%from the finiteness of data translates into uncertainties in parameters.
%In the exposition of the maxent formulation derived from Eq~\ref{eq:Lagrang}, there is no ambiguity in the model parameters because they are fully specified by the constraints calculated from data. In reality, these constraints are typically estimated from a finite sample, and they are noisy.
%The straightforward answer to this problem is to take more data---in a pairwise
%maxent problem, we might insist that we have enough samples to well-constrain
%the correlation between every pair of individuals.  But it is not always possible
%to take enough data.  For instance, in a social system in which we are trying to
%measure stable social structure that lasts on the order of months, there are only
%a finite number of social interactions that occur over those months, which may
%not be enough to tightly constrain parameters. When we fit exactly to constraints measured from a small amount of data, we run the danger of overfitting and poor generalization.
%
%When we find the parameters that minimize the Kullback-Leilber divergence between the model and the data distributions, we are maximizing the likelihood of the data. 
%Sample size influences curvature of the likelihood function around the maximum-likelihood peak.  Specifically, small sample size leads to a smaller curvature, implying a broad peak and larger uncertainty in parameters.
%We can estimate this uncertainty from the data: Assuming that the data is independently and identically distributed, the errors are given by the standard error of the mean of a binomial distribution \mbox{$p_{\rm ij...k} = p(s_{\rm i}=s_{\rm j}=...=s_{\rm k}=1)$} with $R$ data points.
%\begin{align}
%	\delta_{\rm ij...k} &= \sqrt{(1-p_{\rm ij...k})p_{\rm ij...k}/R}.
%\end{align}
%In algorithms that iteratively bring the model closer to the data, we can stop once we have gotten close enough, where close enough is given by the expected fluctuations in the data distribution.
%
%An alternative approach is to regularize the problem as in the regularized mean field or cluster expansion algorithms. Here, we restrict the search space in some principled way so that more complicated
%solutions are disallowed.  We then check that the regularized solutions fit the
%data within expected statistical fluctuations.  If not, a more lax regularization
%can be used to allow more complicated solutions that are able to fit the
%remaining signal in the data. 
%In a Bayesian formulation, this approach is equivalent to including a prior distribution over  parameters that more heavily weights simpler models.
%%\begin{align}
%%    \log p({\rm model}|{\rm data}) &\propto \log p({\rm data}|{\rm model}) + \log p({\rm model})\label{eq:prior}
%%\end{align}
%
%\section*{Model Validation}
%Besides just fitting the model while accounting for the noisiness of the data, how do we know if we need a more complicated model? To answer this question, we need to quantify how well the model fits and the corresponding tradeoff between the complexity and fit of the model.
%
%One way to formalize this is to quantify the tradeoff between a good description of the data with the cost of describing the model. In the Bayesian formulation, we would write this as maximizing the posterior probability of the data given the parameters $\theta$,
%\begin{align}
%    \log p(\theta|\mathcal{D}) &= \log p(\mathcal{D}|\theta) +\log p(\theta) + {\rm const}
%\end{align}
%and the tradeoff manifests as a competition between the likelihood $p(\mathcal{D}|\theta)$ and prior $p(\theta)$. For example, in the regularized mean-field method, the quadratic prior on the couplings favors sparse solutions over more complex solutions that might capture the data better but fail to generalize.
%%in Eq~\ref{eq:prior}. Roughly speaking, we can imagine that specifying the $L$ parameters is localizing a region in a high-dimensional space where each dimension shrinks with the number of data points $K$ as $\sim K^{-1/2}$. The volume of the ``ball'' grows like $K^{-L/2}$, so that the information cost goes like $-\frac{L}{2}\log(K)$ \cite{Lee:2015ev}.
%This picture is formalized by quantities like the Akaike information criterion (AIC) or the Bayesian information criterion (BIC) \cite{Anonymous:mVL3xTtr}.
%
%Although likelihood tells us how well the models are doing relative to one another, it does not tell us in detail how the model is fitting the data.
%Basic checks would be to compare against other correlations in the data.  This could include higher-order correlations than those used to fit the model or other features that are relevant to the question at hand. For instance, an important coarse-grained feature in  macaque conflict is the distribution of fight sizes, which is well reproduced by fitting only pairwise correlations \cite{Daniels:2017cq}.
%% In Figure~??, we show an example of a pairwise maxent model tested against higher order correlations. 
%The most strict check would be to compare the entire model probability distribution with that of the data \cite{Lee:2015ev}, but this is only feasible when the data set is reasonably large.
%
%A summary of how well the model captures the distribution across the entire probability distribution is the multi-information \cite{Nemenman:2004wz,Schneidman:2006he}. If we take a maxent model and add further constraints, the models can be ordered in terms of entropy $S_{\rm 1}>S_{2}>...S_{\rm m}>...>S_{\rm data}$, where the minimum entropy the most constrained model could have is equal to the data. To measure how much correlation in the data our model has captured, we can calculate the multi-information $I_{\rm m} = S_{\rm 1}-S_{\rm m}$, the amount of correlation captured by the model relative to the independent model. The fraction of multi-information captured is $F=I_{\rm m}/I_{\rm data}$. This is a measure of how much of the correlation in the data is captured by the model.\footnote{See Refs \cite{Bialek:2012ueb} and \cite{Lee:2015ev} for details and further references on how to estimate information quantities.}
%
%These tools to quantify model fit then inform the important choice of which constraints to impose. Typically, the approach is to constrain the lowest order interactions that are sufficient to produce collective behavior. The intuition from physics is from the observation that many physical systems are very well described by pairwise interactions \cite{Ranganathan:2007wz}.
%From the model fitting perspective, however, we might choose to constrain other parts of the probability distribution. Ref \cite{Ganmor:2011ct} explores choosing correlations to constrain depending on whether they are signicantly large. Ref \cite{Nemenman:2016kl} discusses how a pairwise model becomes an increasingly effective description of a system with higher order interactions as the system gets larger.


\section*{Implementation and architecture}
The package is divided into three principal modules containing the algorithms for solving the inverse maxent problem ({\tt solvers.py}), the Monte Carlo Markov Chain (MCMC) sampling algorithms ({\tt samplers.py}), and supporting ``utility'' functions for the remaining modules ({\tt utils.py}) as shown in Fig.~\ref{gr:architecture}.
Besides the {\tt utils.py} module, the package is organized around classes that correspond to different algorithms. This class-based structure ensures that the state of the solver or sampler, including the data it was fit to and the current guess for the parameters, are all contained within the instance of the algorithm class. As a result, the current state of work can be saved and moved between workstations using the Python package {\tt dill}.

For the solvers, the different algorithms available are accessible from the {\tt coniii.solvers} module as listed in Fig.~\ref{gr:architecture}. These algorithm classes are all derived from a base {\tt Solver} class as shown in Fig.~\ref{gr:architecture}. In the accompanying release version of the package, v1.0.3, the method {\tt Solver.solve} serves as the interface for solving the inverse maxent problem. To keep the solution algorithms generic enough to solve a variety of different maxent problems, they all require that the user define the maxent model upon instantiation through the definition of keyword arguments like {\tt calc\_observables}. The particular methods required to specify the maxent problem differ by algorithm, but for the pairwise maxent problem we have made it easy by defining those functions as part of the package. These helper functions are available as part of the {\tt utils.py} module and their use is demonstrated in the Jupyter notebook usage guide.

The MCMC sampling algorithms are likewise based on a class architecture derived from {\tt Sampler} as shown in Fig.~\ref{gr:architecture}. Each instance of {\tt Solver} automatically instantiates this class under {\tt Solver.sampler} and wraps calls to it. For the Ising model, this is an instance of {\tt Metropolis}. Other sampling algorithms listed in the {\tt samplers.py} box in Fig.~\ref{gr:architecture} will be released with later versions of this package.

\begin{figure}[tb]
	\includegraphics[width=\linewidth,clip=True,trim=70 450 60 0]{architecture}
	\caption{ConIII architecture. The principal modules are {\tt solvers.py}, {\tt samplers.py}, and {\tt utils.py}. The module {\tt solvers.py} contains classes based on {\tt Solver} that each implement a different algorithm for solving the relevant inverse maxent problem accessible through the method {\tt Solver.solve()}. The {\tt samplers.py} module contains the Metropolis algorithm for Monte Carlo Markov Chain sampling and will support other samplers in future versions (gray font) including Wolff sampling, Swendsen-Wang sampling, and parallel tempering. The {\tt utils.py} module contains supporting functions for the other modules such as the few examples listed. ConIII's modularized structure ensures that contributed algorithms can be appended independently of existing code.}
	\label{gr:architecture}
\end{figure}


\section*{Quality Control}
For checks of basic functionality, the package is released with unit tests that can be run with the Python package \href{https://pytest.org}{pytest}.

The most direct test of the algorithms is to generate a system where the parameters are known, sample from the system to generate a data set, and run the inverse solution to make sure that the correlations and parameters match the known values. With a finite sample, exact correspondence to the correct parameters is not expected although differences should decrease with a larger sample. Furthermore, most of the algorithms only return an approximate solution such that the fidelity of the found parameters to the original ones will depend on the sample size and whether or not the approximation is valid. The Jupyter notebook released with the software provides examples for using the algorithms included in ConIII for a random system of five spins. We recommend that the user run this notebook to check how well different algorithms converge to the solutions depending on the algorithm and sample size.

More importantly, the user can check if the algorithms match the expected correlations closely or not. How one checks the validity of a particular maxent model for data is beyond the scope of this paper, but we point the reader to the appendix of Ref~\cite{Lee:2015ev} where the methodology is explained in detail for a broad audience.

If there are any issues or bugs in the software, we organize improvements and patches through  the \href{https://github.com/eltrompetero/coniii}{GitHub repository} where both issues can be filed and pull requests made. 


\section*{(2) Availability}
\vspace{0.5cm}
\section*{Operating system}
Linux, MacOS, Windows

\section*{Programming language}
Python 3.6, 3.7

\section*{Dependencies}
Python packages multiprocess $\geq$v0.70.5 and $<$v1, numpy, scipy, joblib, matplotlib, numba $\geq$ v0.39.0, dill.

\section*{Software location:}
\begin{description}[noitemsep,topsep=0pt]
	\item[Name:] PyPI
	\item[Persistent identifier:] \url{https://pypi.org/project/coniii/}
	\item[Licence:] MIT License
	\item[Publisher:]  Edward D. Lee
	\item[Version published:] v1.1.1
	\item[Date published:] 12/12/18
\end{description}

{\bf Code repository}

\begin{description}[noitemsep,topsep=0pt]
	\item[Name:] GitHub
	\item[Persistent identifier:] \url{https://github.com/eltrompetero/coniii}
	\item[Licence:] MIT License
	\item[Version published:] v1.1.1
\end{description}

\begin{description}[noitemsep,topsep=0pt]
	\item[Name:] Zenodo
	\item[Persistent identifier:] \url{https://doi.org/10.5281/zenodo.2236632}
	\item[Licence:] MIT License
	\item[Version published:] v1.1.1
\end{description}


\section*{Language}
English


\section*{(3) Reuse potential}
To contribute either an algorithm for the inverse maxent problem or a sampling technique, we suggest following the template for the classes described in the base {\tt Solver} and {\tt Sampler} classes. New algorithms should be filed as a pull request to the \href{https://github.com/eltrompetero/coniii}{GitHub repository} along with an example solution that can be included in the usage guide Jupyter notebook and unit tests.

Documentation for the package is included as part of the GitHub repository and also hosted online at \url{https://eddielee.co/coniii/index.html}.

\section*{Acknowledgements}
We thank the anonymous reviewers for helpful feedback on both the manuscript and the accessibility of the software.


\section*{Funding statement}
EDL was supported by an NSF Graduate Fellowship under grant no. DGE-1650441. This research was supported in part by a congressional research grant provided by the Dirksen Congressional Center.


\section*{Competing interests}
The authors declare that they have no competing interests.

%\section*{Conclusion}
%Maxent provides a coherent framework for modeling a diverse body of biological and social systems.  The approach is structured to consider easily issues of model selection,  embodying a canonical approach for the scientific process: we start with the simplest models and the fewest constraints possible and then increase the complexity of the model until it makes good predictions. In principle, we could add as many constraints as would allow us to better fit the data, but complex models are more computationally expensive and often do not generalize well.  Maxent therefore stops before adding all possible constraints, leaving the model otherwise as structureless as possible.
%In this sense, the maxent approach is a principled method of statistical inference that is relevant beyond its roots in statistical physics.  We have built ConIII as an accessible software package in the hope that it will inspire those unfamiliar with maxent approaches to experiment and apply this technique to their own research questions.

\bibliography{refs}

\vspace{2cm}

\rule{\textwidth}{1pt}


%\section*{How to use}
%Each algorithm class instance must be instantiated with a few parameters that describe the particular maxent problem at hand. For example, all of the classes require a {\tt calc\_observables} keyword argument to be defined. This function will take a 2-dimensional sample (each row is a single state vector and each column is an entry in the state vector) from the ensemble and calculate the observables that correspond to the maxent model of interest. Some algorithms require other functions that are specific both to the problem and the algorithm. For example, the pseudolikelihood algorithm needs the likelihood of a single entry in the state vector ({\it i.e.} a single spin or a single vote) conditional on the rest of the state which depends on the maxent model. We have written helper methods for defining necessary inputs for all the available algorithms for the pairwise maxent problem. Examples of how to use these are given in the accompanying Jupyter notebook.
%
%Once an algorithm class is instantiated, it can be run by calling the {\tt solve} method and by providing the data in the corresponding format for the algorithm. Some algorithms like {\tt MCH} only require the correlations that are constrained in the maxent model. Others like {\tt MPF} require the full data set. 
%
%To check if an algorithm has converged to reasonably close agreement with the data, 
%
%Examples of this workflow for all the provided algorithms are given in the usage guide accompanying the Github package.



%\section*{The microcanonical ensemble and maximum entropy}
%The conventional textbook in statistical mechanics first introduces the concept of entropy as a way of counting the phase volume available to the system at a given energy
%\begin{align}
%	S(E) &= k_B \log\Omega(E)
%\end{align}
%where $k_B$ is Boltzmann's constant and $\Omega$ the number of states between $E$ and $E+\delta E$. Temperature is defined as
%\begin{align}
%	\frac{1}{T} &= \frac{\partial S}{\partial E}
%\end{align}
%
% begins with the concept of a small system coupled to a heat bath. In this limit, we can linearly expand thermodynamic quantities about the energy of the bath $E_{\rm bath} = E-E_s$
%\begin{align}
%	S_{\rm bath}(E-E_s) &\approx S_{\rm bath}(E) -E_s\left.\frac{\partial S}{\partial E_s}\right|_{E-E_s}\\
%		&= S_{\rm bath}(E) -\frac{E_s}{T}\\
%	e^{S_{\rm bath}(E-E_s)/k_B} &= e^{S_{\rm bath}(E)/k_B} e^{-E_s/k_BT}
%\intertext{Since the entropy is proportional to the density of states at a particular energy $E_s$, Eq ?? corresponds to}
%	p(E_s) &= e^{-E_s/k_BT}/Z
%\end{align}
%With partition function $Z$. In other words, the Gibbs measure takes exponential form. Eq ?? can be rewritten in terms of free energy
%\begin{align}
%	\log p(E_s) &= -E_s/k_BT -\log(Z)
%\intertext{Averaging both sides over all energy configurations and rearranging}
%	-k_B T\log Z &= \sum_s p(E_s) E_s +k_BT\sum_s p(E_s)\log p(E_s)
%\intertext{Remembering the fundamental postulate of statistical mechanics that all states are equally likely}
%	-k_BT\log Z &= \br{E_s} - TS\\
%	F &= \br{E_s} -TS
%\end{align}
%This equation tells us how the Helmholtz free energy is related to the internal energy and the entropy of the system.
%
%When we say that free energy is minimized for a system with Hamiltonian $E_s$ at equilibrium, we are equivalently saying that entropy is maximized. Entropy maximization is the crux of statistical physics models.

\end{document}  
