approximate solution as the *Bethe reference point*.

Starting from this point, we want now to add one factor $\psi_{ij}$ to produce the distribution $$\label{eq:1link}
{\mathcal{P}}^{(n+1)}({{\mathbf{x}}}) = {\mathcal{P}}^{(n)}({{\mathbf{x}}})\times \frac{\psi_{ij}(x_i,x_j)}{Z_\psi}$$ with $$Z_\psi = \int dx_idx_j p_{ij}^{(n)}(x_i,x_j)\psi_{ij}(x_i,x_j).$$ The log likelihood corresponding to this new distribution reads $${\mathcal{L}}' = {\mathcal{L}}+ \int d{{\mathbf{x}}}\hat{\mathcal{P}}({{\mathbf{x}}})\log\psi_{ij}(x_i,x_j)-\log Z_\psi.$$ Since the the functional derivative w.r.t. $\psi$ is $$\frac{\partial{\mathcal{L}}'}{\partial \psi_{ij}(x_i,x_j)} = \frac{\hat p(x_i,x_j)}{\psi_{
ij}(x_i,x_j)} - 
\frac{p_{ij}^{(n)}(x_i,x_j)}{Z_\psi},$$ $\forall (x_i,x_j)\in \Omega^2$, the maximum is attained for $$\label{eq:1linkfactor}
\psi_{ij}(x_i,x_j) = \frac{\hat p_{ij}(x_i,x_j)}{p_{ij}^{(n)}(x_i,x_j)}\text{ 
with }  Z_\psi=1,$$ where $p^{(n)}(x_i,x_j)$ is the reference marginal distribution obtained from ${\mathcal{P}}^{(n)}$. The correction to the log-likelihood can then be rewritten as $$\label{eq:onelink}
\Delta {\mathcal{L}}= D_{KL}(\hat p_{ij}\Vert p_{ij}^{(n)}).$$ Sorting all the links w.r.t. this quantity yields the (exact) optimal $1$-link correction to be made. The interpretation is therefore immediate: the best candidate is the one for which the current model yields the joint marginal $p_{ij}^{(n)}$ that is most distant target $\hat p_{ij}$. Note that the update mechanism is indifferent to whether the link has to be added or simply modified.

Iterative proportional scaling for GMRF model selection
-------------------------------------------------------

In the statistics literature this procedure is referred to as the *iterative proportional scaling* (IPS) procedure, originally proposed for contingency table estimations and 