Solving Non-parametric Inverse Problem in Continuous Markov Random Field using Loopy Belief Propagation

TL;DR

This paper presents a novel method for solving the inverse problem in continuous Markov random fields with non-parametric pairwise energy functions, using Bethe approximation and orthonormal function expansion to obtain an analytic solution.

Contribution

It introduces a new approach combining Bethe approximation and orthonormal function expansion to analytically solve the inverse problem in continuous MRFs, which was previously intractable.

Findings

Provides an analytic solution within Bethe approximation framework.

Avoids intractable partition function evaluation using Bethe approximation.

Reduces functional optimization to a manageable function optimization.

Abstract

In this paper, we address the inverse problem, or the statistical machine learning problem, in Markov random fields with a non-parametric pair-wise energy function with continuous variables. The inverse problem is formulated by maximum likelihood estimation. The exact treatment of maximum likelihood estimation is intractable because of two problems: (1) it includes the evaluation of the partition function and (2) it is formulated in the form of functional optimization. We avoid Problem (1) by using Bethe approximation. Bethe approximation is an approximation technique equivalent to the loopy belief propagation. Problem (2) can be solved by using orthonormal function expansion. Orthonormal function expansion can reduce a functional optimization problem to a function optimization problem. Our method can provide an analytic form of the solution of the inverse problem within the framework of Bethe approximation.

Equations

$$\Psi (\mathbf{x}):=-\sum _{i\in V}{\theta }_i(x_i)-\sum _{\{i,j\}\in E}{w}_{\{i,j\}}({\mathbf{x}}_{\{i,j\}}),$$

$$\displaystyle P(\bm{x}):=\frac{1}{Z}\exp\big{(}-\Psi(\bm{x})\big{)},$$

$$\displaystyle Z:=\int_{\mathcal{X}^{n}}\exp\big{(}-\Psi(\bm{x})\big{)}d\bm{x},$$

$$\mathcal{F}[\mathbf{b},\mathbf{\xi }]$$

$$=-\sum _{i\in V}{\int }_{\mathcal{X}}{\theta }_i(x_i)b_i(x_i)d{x}_i-\sum _{\{i,j\}\in E}{\int }_{{\mathcal{X}}^2}{w}_{\{i,j\}}({\mathbf{x}}_{\{i,j\}}){\xi }_{\{i,j\}}({\mathbf{x}}_{\{i,j\}})d{x}_{i}d{x}_j+\sum _{i\in V}(1-|{\partial }_i|){\int }_{\mathcal{X}}{b}_i(x_i)\ln b_i(x_i)d{x}_i$$

$$+\sum _{\{i,j\}\in E}{\int }_{{\mathcal{X}}^2}{\xi }_{\{i,j\}}({\mathbf{x}}_{\{i,j\}})\ln {\xi }_{\{i,j\}}({\mathbf{x}}_{\{i,j\}})d{x}_{i}d{x}_j,$$

$${\int }_{\mathcal{X}}{b}_i(x_i)d{x}_i={\int }_{{\mathcal{X}}^2}{\xi }_{\{i,j\}}({\mathbf{x}}_{\{i,j\}})d{x}_{i}d{x}_j=1,$$

$${\int }_{\mathcal{X}}{\xi }_{\{i,j\}}({\mathbf{x}}_{\{i,j\}})d{x}_i$$

$${\int }_{\mathcal{X}}{\xi }_{\{i,j\}}({\mathbf{x}}_{\{i,j\}})d{x}_j$$

$$m_{i\to j}(x_j)=\frac{1}{Z_{i\to j}}{\int }_{\mathcal{X}}{\pi }_{i\setminus j}(x_i)e^{w_{\{i,j\}}({\mathbf{x}}_{\{i,j\}})}d{x}_i,$$

$$Z_{i\to j}:={\int }_{{\mathcal{X}}^2}{\pi }_{i\setminus j}(x_i)e^{w_{\{i,j\}}({\mathbf{x}}_{\{i,j\}})}d{x}_{i}d{x}_j$$

$${\pi }_{i\setminus j}(x_i):=\frac{e^{{\theta }_i(x_i)}{\prod }_{k\in {\partial }_i\setminus \{j\}}{m}_{k\to i}(x_i)}{{\int }_{\mathcal{X}}{e}^{{\theta }_i(x_i)}{\prod }_{k\in {\partial }_i\setminus \{j\}}{m}_{k\to i}(x_i)d{x}_i}.$$

$$b_i(x_i)$$

$${\xi }_{\{i,j\}}({\mathbf{x}}_{\{i,j\}})$$

$$l[\mathbf{\theta },\mathbf{w}]:=\frac{1}{N}\sum _{\mu =1}^N\ln P({\mathbf{x}}^{(\mu )}),$$

$$l[\mathbf{\theta },\mathbf{w}]$$

$$l_{\mathrm{Bethe}}[\mathbf{\theta },\mathbf{w}]:=-\frac{1}{N}\sum _{\mu =1}^N\Psi ({\mathbf{x}}^{(\mu )})+\underset{\mathbf{b},\mathbf{\xi }}{\min }\mathcal{F}[\mathbf{b},\mathbf{\xi }].$$

$$\frac{\delta l_{\mathrm{Bethe}}[\mathbf{\theta },\mathbf{w}]}{\delta {\theta }_i(x_i)}$$

$$\frac{\delta l_{\mathrm{Bethe}}[\mathbf{\theta },\mathbf{w}]}{\delta w_{\{i,j\}}({\mathbf{x}}_{\{i,j\}})}$$

$$-{\xi }_{\{i,j\}}({\mathbf{x}}_{\{i,j\}}),$$

$$\frac{1}{N}\sum _{\mu =1}^{N}f({\mathrm{x}}_i^{(\mu )})={\int }_{\mathcal{X}}f(x_i)b_i(x_i)d{x}_i$$

$$\frac{1}{N}\sum _{\mu =1}^{N}g({\mathrm{x}}_i^{(\mu )},{\mathrm{x}}_j^{(\mu )})={\int }_{\mathcal{X}}g(x_i,x_j){\xi }_{\{i,j\}}({\mathbf{x}}_{\{i,j\}})d{x}_{i}d{x}_j$$

$${\int }_{\mathcal{X}}{\varphi }_s(x){\varphi }_t(x)dx={\delta }_{s,t},$$

$$f(x)=\sum _{s=0}^{\infty }{\alpha }_s{\varphi }_s(x),$$

$${\alpha }_s={\int }_{\mathcal{X}}f(x){\varphi }_s(x)dx.$$

$${\varphi }_0(x)=\frac{1}{\sqrt{\chi }},\chi :=\beta -\alpha .$$

$${\int }_{\mathcal{X}}{\varphi }_s(x)dx=\{\begin{array}{ll}\sqrt{\chi }& s=0\\ 0& s>0\end{array}.$$

$${\theta }_i(x_i)$$

$$w_{\{i,j\}}({\mathbf{x}}_{\{i,j\}})$$

$$\displaystyle=\frac{1}{\sqrt{\chi}}\sum_{s=1}^{\infty}\big{(}J_{\{i,j\}}^{(s,0)}\phi_{s}(x_{i})+J_{\{i,j\}}^{(0,s)}\phi_{s}(x_{j})\big{)}$$

$$+\sum _{s,t=1}^{\infty }{J}_{\{i,j\}}^{(s,t)}{\varphi }_s(x_i){\varphi }_t(x_j)+\mathrm{constant},$$

$$h_i^{(s)}$$

$$J_{\{i,j\}}^{(s,t)}$$

$$\Psi (\mathbf{x})={\Psi }^{†}(\mathbf{x};\mathbf{H},\mathbf{J})+C_0,$$

$${\Psi }^{†}(\mathbf{x};\mathbf{H},\mathbf{J})$$

Full text

Markov random field, Boltzmann machine, continuous variable, loopy belief propagation, nonparametric estimation, orthonormal function expansion

††preprint: APS/123-QED

I Introduction

Boltzmann machine learning, which is known as the inverse Ising problem in statistical mechanics, is one of the important problems in the statistical machine learning field and has a long history. Suppose that we have sample points, i.e., data points, stochastically generated from an unknown distribution (referred to as a generative model). The task of statistical machine learning is to specify the unknown distribution using only the sample points. In standard Boltzmann machine learning, we assume that the generative model that generates data points is an Ising model, and prepare an Ising model (referred to as the learning model) with controllable parameters, e.g., external fields and exchange interactions. The Boltzmann machine learning is achieved by optimizing the values of the controllable parameters in the learning model through maximum likelihood estimation.

Unfortunately, we cannot perform Boltzmann machine learning exactly because of the computational cost. Therefore, many approximations for Boltzmann machine learning have been proposed. In particular, approximations based on mean-field methods have been developed in the field of statistical mechanics Roudi et al. (2009): mean-field approximation Kappen and Rodríguez (1998), Bethe approximation Parise and Welling (2009); Yasuda and Horiguchi (2006); Yasuda and Tanaka (2009); Mézard and Mora (2009); Marinari and Kerrebroeck (2010); Ricci-Tersenghi (2012); Nguyen and Berg (2012); Furtlehner (2013), Plefka expansion Tanaka (1998); Sessak and Monasson (2009), and so on. In many of these methods, we can obtain the solution to the maximum likelihood estimation analytically. However, they are applicable to only an Ising-type learning model, that is, the variables in the model are binary and the energy function of the model is a quadratic form of the variables.

We proposed a method for a more general situation that uses Bethe approximation and orthonormal function expansion Yasuda et al. (2012). Using the method, we can solve the inverse problem with general pair-wise Markov random fields and obtain the solution analytically. However, this method cannot be applied to Markov random fields with continuous variables.

In this paper, we propose a method for solving the inverse problem in general pair-wise Markov random fields with continuous variables, which is an extension of our previous method Yasuda et al. (2012). The proposed method can give us the analytical solution of the inverse problem. This is the main contribution of this paper. In this paper, we refer to a pair-wise Markov random field with continuous variables as a continuous Markov random field (CMRF).

The remainder of this paper is organized as follows. In Sec. II, we explain loopy belief propagation (LBP) in a CMRF. LBP is equivalent to Bethe approximation Yedidia et al. (2005); Pelizzola (2005). We formulate the inverse problem in a CMRF in Sec. III, as well as its Bethe approximation. Our method is shown in Sec. IV. In this section, we derive the solution to the inverse problem using the Bethe approximation shown in Sec. III. Since the solution is obtained in the form of infinite series, it cannot be implemented as it is. We describe a means of implementing our method and show the results of numerical experiments in Sec. V. We conclude the paper with some remarks in Sec. VI.

II Formalism of Loopy Belief Propagation in Continuous Markov Random Field

Consider an undirected graph $G(V,E)$, where $V=\{1,2,\ldots,n\}$ is the set of nodes and $E$ is the set of undirected links. We denote the link between nodes $i$ and $j$ by $\{i,j\}$. Because the links have no direction, $\{i,j\}$ and $\{j,i\}$ indicate the same link. On the undirected graph, we define the non-parametrized pair-wise energy function as

[TABLE]

where $\theta_{i}(x_{i})$ is the energy on node $i$, $w_{\{i,j\}}(\bm{x}_{\{i,j\}})$ is the energy on link $\{i,j\}$, and $\bm{x}_{\{i,j\}}=\bm{x}_{\{j,i\}}=\{x_{i},x_{j}\}$. We regard $w_{\{i,j\}}(\bm{x}_{\{i,j\}})$ as the same function as $w_{\{j,i\}}(\bm{x}_{\{i,j\}})$. With the energy function, we define the CMRF as

[TABLE]

where $\bm{x}=\{x_{i}\in\mathcal{X}\mid i\in V\}$ represents the continuous random variables over the continuous space $\mathcal{X}$, and $Z$ is the partition function defined as

[TABLE]

where $\int_{\mathcal{X}^{n}}f(\bm{x})d\bm{x}$ denotes the multiple integration over whole variables, $\int_{\mathcal{X}}\cdots\int_{\mathcal{X}}f(\bm{x})dx_{1}\cdots dx_{n}$, and $\int_{\mathcal{X}}$ denotes the integral over $\mathcal{X}$. $\theta_{i}(x_{i})$ and $w_{\{i,j\}}(\bm{x}_{\{i,j\}})$ are arbitrary functions of the assigned variables.

Given the CMRF, it is difficult to evaluate its marginal distributions because of the existence of intractable multiple integration. LBP is one of the most effective methods for approximately evaluating marginal distributions and is the same as Bethe approximation in statistical mechanics. LBP can be obtained from the minimum condition of the variational Bethe free energy of the CMRF in Eq. (2). We denote the marginal distribution over $x_{i}$ by $b_{i}(x_{i})$ and that over $x_{i}$ and $x_{j}$, which are neighboring pair of nodes, by $\xi_{\{i,j\}}(\bm{x}_{\{i,j\}})$. These marginal distributions are sometimes called beliefs in the context of LBP. We regard $\xi_{\{j,i\}}(\bm{x}_{\{i,j\}})$ as the same belief as $\xi_{\{i,j\}}(\bm{x}_{\{i,j\}})$. In the context of the cluster variation method Kikuchi (1951); Yedidia et al. (2005), the variational Bethe free energy of the CMRF is expressed as

[TABLE]

where $\partial_{i}=\{j\mid\{i,j\}\in E\}$ is the set of nodes connected to node $i$. The variational Bethe free energy is regarded as the functional with respect to $\bm{b}=\{b_{i}(x_{i})\mid i\in V\}$ and $\bm{\xi}=\{\xi_{\{i,j\}}(\bm{x}_{\{i,j\}})\mid\{i,j\}\in E\}$. The beliefs, that minimize the variational Bethe free energy, are regarded as the Bethe approximation of the corresponding marginal distributions. From the extremal condition of the variational Bethe free energy under the normalizing constraints,

[TABLE]

and the marginalizing constraints,

[TABLE]

we obtain the message-passing equation (MPE)

[TABLE]

where the constant $Z_{i\to j}$ is frequently set to

[TABLE]

to normalize the messages. The distribution $\pi_{i\setminus j}(x_{i})$ is defined as

[TABLE]

The quantity $m_{i\to j}(x_{j})$ is the normalized message (or the effective field) from node $i$ to node $j$, which is non-negative and originates from the Lagrange multipliers appearing in the conditional minimization of the variational Bethe free energy. The two different messages, $m_{i\to j}(x_{j})$ and $m_{j\to i}(x_{i})$, are defined on link $\{i,j\}$. The beliefs (the approximate marginal distributions) are computed from the messages as

[TABLE]

In principle, by solving the MPE in Eq. (7), we can compute the one- and two-variable marginal distributions using Eqs. (10) and (11). However, finding the functional forms of the messages is not straightforward, because the messages are continuous functions over $\mathcal{X}$, and therefore, the MPE we have to solve is an integral equation. Some methods that are based mainly on a stochastic method have been developed for approximately solving the MPE Sudderth et al. (2003); Ihler and McAllester (2009); Noorshams and Wainwright (2013).

III Inverse Problem in Continuous Markov Random Field

In this section, we consider the inverse problem, in other words, the machine learning problem, for the CMRF in Eq. (2). The inverse problem for the CMRF can be solved by maximum likelihood estimation. Given $N$ data points $\mathcal{D}=\{\mathbf{x}^{(\mu)}\in\mathcal{X}^{n}\mid i=1,2,\ldots,N\}$, we define the log-likelihood functional as

[TABLE]

where $\bm{\theta}$ and $\bm{w}$ are the set of functions $\theta_{i}(x_{i})$ and $w_{\{i,j\}}(\bm{x}_{\{i,j\}})$ respectively in the exponent in Eq. (2). The goal of the maximum likelihood estimation is to find the functions $\bm{\theta}$ and $\bm{w}$ that maximize the log-likelihood functional. Eq. (12) can be rewritten as

[TABLE]

However, the maximization problem of the log-likelihood functional is intractable because of the existence of the partition function.

To avoid evaluating the intractable partition function, we approximate the log-likelihood functional using LBP, i.e., Bethe approximation. The Bethe approximation of the log-likelihood functional in Eq. (13) can be expressed by using the variational Bethe free energy shown in Eq. (3) as

[TABLE]

We refer to this as the Bethe log-likelihood functional. The main purpose of this study was to maximize the Bethe log-likelihood functional with respect to the functions $\bm{\theta}$ and $\bm{w}$. The solution obtained by maximizing Eq. (14), of course coincides to that obtained by the true maximum likelihood estimation when the CMRF has a tree structure, because Bethe approximation is exact in tree systems. However, the maximization of the Bethe log-likelihood functional is not straightforward for the following reasons. The variations of the functional with respect to $\bm{\theta}$ and $\bm{w}$ are

[TABLE]

where $b_{i}(x_{i})$ and $\xi_{\{i,j\}}(\bm{x}_{\{i,j\}})$ are the beliefs minimizing the variational Bethe free energy, in other words, the solution to the LBP presented in the previous section. This variation means that we have to find $\bm{\theta}$ and $\bm{w}$ that satisfy the relations

[TABLE]

and

[TABLE]

for any test functions $f(x_{i})$ and $g(x_{i},x_{j})$. Thus, if we could obtain the solution of the LBP, by using a method that has already proposed Sudderth et al. (2003); Ihler and McAllester (2009); Noorshams and Wainwright (2013), the solution to the maximization of the Bethe log-likelihood functional is not immediately obtained.

IV Proposed Method

In this section, we propose a method to solve the maximization problem of the Bethe log-likelihood function in Eq. (14) in terms of orthonormal function expansion. Via orthonormal function expansion, we can reduce the functional maximization problem in the previous section to a tractable function maximization problem. The basic idea of our method is similar to that presented in our previous paper Yasuda et al. (2012).

IV.1 Orthonormal Function System

Before deriving our method, we introduce an orthonormal function system $\{\phi_{s}(x)\mid s=0,1,2,\ldots\}$ over $\mathcal{X}$ satisfying

[TABLE]

where $\delta_{s,t}$ is the Kronecker delta function. By using the orthonormal function system, function $f(x)$ over $\mathcal{X}$ is expanded as

[TABLE]

where the expanding coefficients are given by

[TABLE]

The orthonormal function expansion in Eq. (18) plays an important role in our method.

In the following, we assume that $\mathcal{X}$ is the finite space, $\mathcal{X}=[\alpha,\beta]$, and that $\phi_{0}(x)$ is constant over $\mathcal{X}$, i.e.,

[TABLE]

From Eqs. (17) and (20), we have

[TABLE]

Examples of this orthonormal function are described in Appendix A. We use Eqs. (17), (20), and (21) frequently throughout the paper.

The orthonormal function expansion introduced in this section plays a central role in our proposed method described in the following. However, a similar idea can be useful for solving the LBP in Sec. II. Indeed, a method for solving the LBP was proposed by using orthonormal function expansion Noorshams and Wainwright (2013).

IV.2 Variational Bethe Free Energy with Orthonormal Function Expansion

First, we rewrite the CMRF in Eq. (2) by expanding $\bm{\theta}$ and $\bm{w}$. By using the orthonormal function expansion in Eq. (18), the functions $\bm{\theta}$ and $\bm{w}$ can be expanded as follows.

[TABLE]

and

[TABLE]

where, from Eq. (19), the expanding coefficients are

[TABLE]

It is noteworthy that, from the symmetric property of $w_{\{i,j\}}(\bm{x}_{\{i,j\}})$, $J_{\{i,j\}}^{(s,t)}=J_{\{j,i\}}^{(t,s)}$ is satisfied. In Eqs. (22) and (23), Eq. (20) is used. Using Eqs. (22) and (23), we can rewrite the energy function in Eq. (1) as

[TABLE]

where

[TABLE]

and

[TABLE]

The constant $C_{0}$ in Eq. (26) originates from the constants in Eqs. (22) and (23). Therefore, using the new energy function, the CMRF in Eq. (2) can be rewritten as

[TABLE]

This rewriting makes the CMRF the parametric model, parameterized by $\bm{H}=\{H_{i}^{(s)}\mid i\in V,s\geq 1\}$ and $\bm{J}=\{J_{\{i,j\}}^{(s,t)}\mid\{i,j\}\in E,s\geq 1,t\geq 1\}$. In Eq. (28), the constant in Eq. (26) is neglected, because it is irrelevant to the distribution.

Now, we introduce the orthonormal function expansions of the beliefs in the variational Bethe free energy, as follows.

[TABLE]

From Eq. (19), the expanding coefficients are

[TABLE]

The beliefs must satisfy the normalizing constraints in Eq. (4) and the marginalizing constraints in Eqs. (5) and (6). From Eqs. (4), (20), (31), and (32), we have

[TABLE]

From Eqs. (5), (6), (20), (31), and (32), we obtain

[TABLE]

From the above equations, the beliefs in Eqs. (29) and (30) can be expressed as

[TABLE]

where $\bm{c}=\{c_{i}^{(s)}\mid i\in V,s\geq 1\}$ and $\bm{d}=\{d_{\{i,j\}}^{(s,t)}\mid\{i,j\}\in E,s\geq 1,t\geq 1\}$. By using Eq. (21), one can confirm that the beliefs in Eqs. (33) and (34) satisfy the normalization constraints and the marginal constraints for any $\bm{c}$ and $\bm{d}$.

From Eqs. (27), (33), and (34), in the same way as in Eq. (3), we formulate the variational Bethe free energy for the CMRF in Eq. (28) as

[TABLE]

For specific $\bm{\theta}$ and $\bm{w}$, this variational Bethe free energy is not the functional, but the function of $\bm{c}$ and $\bm{d}$. The variational Bethe free energy in Eq. (35) coincides with that in Eq. (3), except for the irrelevant constant neglected in Eq. (28), i.e.,

[TABLE]

As mentioned above, the beliefs in Eqs. (33) and (34) satisfy the normalization constraints and the marginal constraints for any $\bm{c}$ and $\bm{d}$, so that we can minimize $\mathcal{F}(\bm{c},\bm{d})$ with no constraint. At the minimum point of $\mathcal{F}(\bm{c},\bm{d})$, $\bm{c}$ and $\bm{d}$ satisfy

[TABLE]

Eqs. (37) and (38) are derived from the extremal condition of Eq. (35) with respect to $c_{i}^{(s)}$ and $d_{\{i,j\}}^{(s,t)}$, respectively. In the derivation of these equations, we used Eq. (21).

IV.3 Maximization of the Bethe Log-likelihood Function

By using the new energy function in Eqs. (27) and the variational Bethe free energy in Eq. (35), the Bethe log-likelihood functional in Eq. (14) is represented as

[TABLE]

This is the function with respect to $\bm{H}$ and $\bm{J}$ and we refer to this function as the Bethe log-likelihood function. Thus, the functional optimization problem of the maximum likelihood estimation is reduced to the function optimization problem. The Bethe log-likelihood function is equivalent to the Bethe log-likelihood functional in Eq. (14), because, from Eqs. (26) and (36),

[TABLE]

Therefore, the maximization of the the Bethe log-likelihood function with respect to $\bm{H}$ and $\bm{J}$ is equivalent to the maximization of the Bethe log-likelihood functional with respect to $\bm{\theta}$ and $\bm{w}$. At the maximum point of the Bethe log-likelihood function, we have equations for the expanding coefficients in Eqs. (33) and (34) as

[TABLE]

where $\langle f(\bm{x})\rangle_{\mathcal{D}}:=N^{-1}\sum_{\mu=1}^{N}f(\mathbf{x}^{(\mu)})$ is the sample average over data points $\mathcal{D}$. Coefficients $\hat{c}_{i}^{(s)}$ and $\hat{d}_{\{i,j\}}^{(s,t)}$ are the solutions to the minimization of the variational Bethe free energy in Eq. (35), that is, the solutions to Eqs. (37) and (38). In the following, we denote the beliefs, the coefficients of which are fixed by Eqs. (40) and (41), by $\hat{\bm{b}}$ and $\hat{\bm{\xi}}$, i.e., $\hat{b}_{i}(x_{i}):=b_{i}(x_{i}\mid\hat{\bm{c}})$ and $\hat{\xi}_{\{i,j\}}(\bm{x}_{\{i,j\}}):=\xi_{\{i,j\}}(\bm{x}_{\{i,j\}}\mid\hat{\bm{c}},\hat{\bm{d}})$.

By substituting Eqs. (40) and (41) into Eqs. (37) and (38), we can obtain the solution, $\hat{\bm{H}}$ and $\hat{\bm{J}}$, to the maximization of the Bethe log-likelihood function in Eq. (39), and then, identify the energy function $\Psi^{\dagger}(\bm{x}\mid\hat{\bm{H}},\hat{\bm{J}})$. It should be noted that the solution obtained by our method satisfies Eqs. (15) and (16), which is easily confirmed as follows. A test function $f(x)$ is expanded as in Eq. (18). Therefore, the left side of Eq. (15) is

[TABLE]

On the other side, the right hand side of Eq. (15) is

[TABLE]

where we use Eq. (33). From these equations and Eq. (40), the solution obtained by our method satisfying Eq. (15) is confirmed. Similarly, we can verify the equality in Eq. (16).

By using the method described above, within the framework of Bethe approximation we can identify the functional form of the energy function through the use of the given $N$ data points, and then obtain the resulting CMRF as

[TABLE]

Unfortunately, one cannot computationally treat the infinite series in Eqs. (37) and (38). Thus, in practice, we truncate the infinite series and approximate them by a finite series obtained by the truncation. The details of this approximation are described in Sec. V.1.

The proposed method includes the integration procedures (cf. Eqs. (37) and (38)). The following rewriting allows us to identify the functional form of $\Psi^{\dagger}(\bm{x}\mid\hat{\bm{H}},\hat{\bm{J}})$ without the integration procedures. We now consider the energy function defined by

[TABLE]

This energy function satisfies the relation

[TABLE]

where $C_{1}$ is the constant unrelated to $\bm{x}$ (see Appendix B). From the relation and Eq. (42), we obtain the CMRF determined by the Bethe approximation of the MLE as

[TABLE]

and obtain the energy function in the form of Eq. (43).

V Implementation

V.1 Approximation for Implementation

The beliefs, $\hat{b}_{i}(x_{i})$ and $\hat{\xi}_{ij}(x_{i},x_{j})$, are expressed by an infinite series, as shown in Eqs. (33) and (34). Because an infinite series is not implementable, we approximate them by the truncation up to a finite order:

[TABLE]

and

[TABLE]

where the positive integer $K$ controls the order of the approximation. In the limit of $K\to\infty$, $\hat{b}_{i}^{(K)}(x_{i})$ and $\hat{\xi}_{\{i,j\}}^{(K)}(\bm{x}_{\{i,j\}})$ coincide with $\hat{b}_{i}(x_{i})$ and $\hat{\xi}_{\{i,j\}}(\bm{x}_{\{i,j\}})$, respectively. The approximate beliefs in Eqs. (46) and (47) are normalized for any $K>0$.

Because of the above truncating approximation, the non-negativity of the beliefs may not be retained. Thus, to preserve the positivity of the beliefs, we have to make a further approximation to them. For a small positive value $\varepsilon$, we define distributions

[TABLE]

and regard the cut-off distribution as the approximation of $\hat{b}_{i}^{(K)}(x_{i})$. If $\hat{b}_{i}^{(K)}(x_{i})\geq\varepsilon$ over $\mathcal{X}$, $\tilde{b}_{i}^{(K)}(x_{i})=\hat{b}_{i}^{(K)}(x_{i})$. In a similar manner, we approximate $\hat{\xi}_{\{i,j\}}^{(K)}(\bm{x}_{\{i,j\}})$ by

[TABLE]

By using Eqs. (48) and (49) instead of $\hat{b}_{i}(x_{i})$ and $\hat{\xi}_{ij}(x_{i},x_{j})$, the CMRF in Eq. (45) is approximated by

[TABLE]

where

[TABLE]

is the approximation of Eq. (43). Constant $C_{2}$, which originates from the denominators of Eqs. (48) and (49), is negligible in Eq. (50).

The procedure of our method is summarized as follows. First, given $\mathcal{D}$ we compute $\hat{\bm{c}}$ and $\hat{\bm{d}}$ in Eqs. (40) and (41). Then, using $\hat{\bm{c}}$ and $\hat{\bm{d}}$, we compute $\tilde{b}_{i}^{(K)}(x_{i})$ and $\tilde{\xi}_{\{i,j\}}^{(K)}(\bm{x}_{\{i,j\}})$ in Eqs. (48) and (49), and then $\tilde{\Psi}_{K}^{\ddagger}(\bm{x})$ in Eq. (51) for certain $K>0$ and $\varepsilon>0$. Finally, we obtain the CMRF determined in our method by Eq. (50), and regard the CMRF as the solution to the inverse problem.

If one wants to obtain coefficients $\hat{\bm{H}}$ and $\hat{\bm{J}}$, one can approximately obtain them by using Eqs. (37) and (38) together with $\tilde{b}_{i}^{(K)}(x_{i})$ and $\tilde{\xi}_{\{i,j\}}^{(K)}(\bm{x}_{\{i,j\}})$ instead of $\hat{b}_{i}(x_{i})$ and $\hat{\xi}_{\{i,j\}}(\bm{x}_{\{i,j\}})$.

V.2 Numerical Experiment

Let us consider CMRF $P_{\mathrm{gen}}(\bm{x})\propto\exp(-\Psi_{\mathrm{gen}}(\bm{x}))$ on an undirected graph $G(V,E)$ with $n=9$, where the energy function is defined as

[TABLE]

$\mu=(\beta-\alpha)/2$, and $\mathcal{X}=[\alpha,\beta]=[0,1]$. Suppose that the CMRF is the generative model lying behind the data points in our numerical experiments. We generate $N$ data points, $\mathcal{D}$, from the generative model by the Markov chain Monte Carlo method, and then solve the inverse problem by the method proposed in the previous section using $\mathcal{D}$. In the following experiments, we supposed that the CMRF used in solving the inverse problem has the same graph structure as the generative CMRF and we used Eq. (60) as the orthonormal function system in our method.

In the first experiment, we supposed that the generative CMRF is defined on a 1D chain graph. Because Bethe approximation gives exact solutions in systems with no loops, our method described in Sec. IV.3 provides the true solution to the maximum likelihood estimation for the true log-likelihood in Eq. (12). Given $\mathcal{D}$, we computed $\tilde{P}_{K}(\bm{x})$ in Eq. (50) for a certain $K>0$ and $\varepsilon=0.0001$ by following the procedure described in Sec. V.1. Fig. 1 shows the log-likelihood defined by

[TABLE]

against various $K$. In the computation of $\tilde{l}_{K}$, we approximately evaluated the partition function in $\tilde{P}_{K}(\bm{x})$ by the Monte Carlo integration:

[TABLE]

where $\mathbf{y}^{(m)}$ is the $m$-th sampled point drawn from the unique distribution over $\mathcal{X}$ and $M=20000$. Note that let $\tilde{P}_{K}(\mathbf{x}^{(\mu)})$ be the unique distribution, $\tilde{P}_{K}(\mathbf{x}^{(\mu)})=1$, when $K=0$ in this experiment. The log-likelihood represents the fitness of the solution, $\tilde{P}_{K}(\bm{x})$, to $\mathcal{D}$. A solution that gives a higher value of $\tilde{l}_{K}$ fits the data set better.

$\tilde{l}_{K}$ increases with the increase in the value of $K$, as shown in Fig. 1. This is because a larger value of $K$ increases the number of controllable parameters and increases the flexibility of the model. A more flexible model fits the data set better. In the plot in Fig. 1 and the following plots, since the error bars (the standard deviations) are too small to be visible, we do not show them.

Is the solution with larger $K$ always better? The answer is no in general. The important purpose of the inverse problem is to reconstruct the generative model using the given data set. It is known that an over-fit to the data set frequently degrades the quality of the reconstruction, because a finite size data set includes noise. We measure the quality of the reconstruction, referred to as the generalization error, by the Kullback-Leibler divergence (KLD) defined as

[TABLE]

The solution that gives a smaller value of $\mathcal{K}_{K}$ constitutes a better reconstruction of the generative model.

Fig. 2 shows the KLD for various values of $K$. The KLDs are approximately evaluated by a certain Monte Carlo integration method. In the perspective of the generalization error, the optimal value of $K$ is $K=2$ when $N=1000$ and $K=4$ when $N=10000$.

In practice, we cannot compute $\mathcal{K}_{K}$, because the generative model is unknown. The Akaike information criterion (AIC) is one of the most useful criteria of the generalization error Akaike (1973). The AIC is defined as

[TABLE]

where $\tilde{l}_{K}$ is the log-likelihood defined in Eq. (53) and

[TABLE]

is the number of controllable parameters. In the context of the AIC, the model that minimizes the AIC is the best in the perspective of the generalization error. Fig. 3 shows the AIC for various $K$.

We confirm that the AIC is minimized at $K=2$ when $N=1000$ and at $K=4$ when $N=10000$ and that these are consistent with the results in Fig. 2.

In the next experiment, we supposed the generative CMRF, $P_{\mathrm{gen}}(\bm{x})$, is defined on a $3\times 3$ square grid graph, and performed the same numerical experiments as those described above. The log-likelihood, the KLD, and the AIC are shown in Figs. 4–6, respectively. We observe the results similar to those of the first experiment.

VI Conclusion

In this paper, we proposed a method for the inverse problem in the CMRF with the non-parametrized pair-wise energy function shown in Eq. (1) that uses LBP and orthonormal function expansion. As shown in Sec. IV.3, our method can provide the analytic solution to the inverse problem in the form of an infinite series. Since one cannot treat the infinite series computationally, we proposed further approximations, the truncation approximation in Eqs. (46) and (47) and the cut-off approximation in Eqs. (48) and (49), for the implementation of our method, in Sec. V.1. The numerical results for artificial data were shown in Sec. V.2. From the numerical results, we observed that the optimal value of truncation order $K$ could be found by the AIC.

However, our method still has a strong limitation, that is, the sample space of the variables is a finite space, $\mathcal{X}=[\alpha,\beta]$. This limitation was required to impose the normalization constraints and the marginal constraints on the beliefs in Eqs. (29) and (30) for any $\bm{c}$ and $\bm{d}$, and this property is quite important for our derivation. Thus, the extension to the case where $\mathcal{X}$ is an infinite space may require an approach different from that used in the current study. The extension of our method to such a case will be addressed in our future works.

Appendix A Examples of Orthonormal Function System

In Sec. IV.1, we assume $\mathcal{X}=[\alpha,\beta]$ and $\phi_{0}(x)=1/\sqrt{\chi}$. One possible choice of the orthonormal function set is the normalized Legendre polynomial, $\phi_{s}(x)=L_{s}(x)$, defined as

[TABLE]

The normalized Legendre polynomial satisfies Eq. (17) on $\mathcal{X}=[-1,1]$ and it satisfies the recursion formula

[TABLE]

for $s\geq 1$, where and $L_{0}(x)=1/\sqrt{2}$ and $L_{1}(x)=\sqrt{3/2}x$. Another possible choice is $\phi_{s}(x)=C_{s}(x)$, where

[TABLE]

This is the orthonormal function system on $\mathcal{X}=[0,\pi]$.

It is noteworthy that, by using linear transformation, we can obtain a new orthonormal function system $\{\tilde{\phi}_{s}(x)\}$ on $\mathcal{X}=[\gamma,\delta]$ from an orthonormal function system $\{\phi_{s}(x)\}$ on $\mathcal{X}=[\alpha,\beta]$ as

[TABLE]

where $\tilde{\chi}:=\delta-\gamma$. The new function system $\{\tilde{\phi}_{s}(x)\}$ satisfies Eq. (17) on $\mathcal{X}=[\gamma,\delta]$, and $\tilde{\phi}_{0}(x)=1/\sqrt{\tilde{\chi}}$, when $\phi_{0}(x)=1/\sqrt{\chi}$.

Appendix B The Relation in Eq. (44)

By using the same technique as in to Eqs. (22) and (23), $\ln\hat{b}_{i}(x_{i})$ and $\ln\hat{\xi}_{\{i,j\}}(\bm{x}_{\{i,j\}})$ are expanded as

[TABLE]

and

[TABLE]

respectively, where $A_{i}^{(s)}$ and $B_{\{i,j\}}^{(st)}$ are the expanding coefficients defined by

[TABLE]

From Eqs. (37), (38), and $B_{\{i,j\}}^{(s,t)}=B_{\{j,i\}}^{(t,s)}$, we have

[TABLE]

From Eqs. (62)–(65), the orthonormal function expansion of the the energy function in Eq. (43) is written as

[TABLE]

where $C_{1}$ is the constant originates from the constants in Eqs. (62) and (63). From this equation we obtain Eq. (44).

Acknowledgment

This work was partially supported by JST CREST Grant Number JPMJCR1402 and by JSPS KAKENHI Grant Numbers 15K00330, 15H03699, and 15K20870.