MLE-induced Likelihood for Markov Random Fields
Jie Liu Hao Zheng


The effects of kinetic instabilities on the electron cyclotron emission from runaway electrons
    Max E. Austin
    December 30, 2023
==============================================================================================




Due to the intractable partition function, the exact likelihood function for a Markov random field (MRF), in many situations, can only be approximated.
Major approximation approaches include pseudolikelihood <cit.> and Laplace approximation <cit.>.
In this paper, we propose a novel way of approximating the likelihood function through first approximating the marginal likelihood functions of individual parameters and then reconstructing the joint likelihood function from these marginal likelihood functions.
For approximating the marginal likelihood functions, we derive a particular likelihood function from a modified scenario of coin tossing which is useful for capturing how one parameter interacts with the remaining parameters in the likelihood function.
For reconstructing the joint likelihood function, we use an appropriate copula to link up these marginal likelihood functions.
Numerical investigation suggests the superior performance of our approach.
Especially as the size of the MRF increases, both the numerical performance and the computational cost of our approach remain consistently satisfactory, whereas Laplace approximation deteriorates and pseudolikelihood becomes computationally unbearable.




§ INTRODUCTION


Suppose that we observe n i.i.d. data points 𝕏={ x^1,..., x^n} sampled from a Markov random field (MRF) on X∈𝒳^d (𝒳 is a discrete space) which is parameterized by θ. 
The likelihood function L(θ|𝕏) plays a key role in statistical inference of an MRF, including parameter estimation, model comparison, simulation practice and others. 
However due to the intractable partition function, the exact likelihood function, in many situations, can only be approximated. 
Widely used approximation approaches include pseudolikelihood <cit.>, Laplace approximation <cit.> and their variations, but these methods may suffer from either precision loss or a heavy computational burden, especially when we are dealing with large-scale MRFs.

Aiming at a fast and reliable approximation of the MRF likelihood function, we propose a two-step method as follows. 
First, we focus on the individual dimensions of the parameter space, and approximate the marginal likelihood <cit.> of the individual parameters. 
We derive a particular likelihood function from a modified scenario of coin tossing, and use this function for marginal likelihood approximation. 
In the second step, we reconstruct the joint likelihood function by linking up these marginal likelihood functions. 
We refer to a popular technique in high-dimensional statistical inference — copula — to reconstruct the joint likelihood. 
Because the first step of our approximation depends on the availability of the maximum likelihood estimate (MLE) of the parameters, we term the achieved likelihood function MLE-induced likelihood.

The rationale of our two-step approximation of MRF likelihood is as follows. 
In the first step of approximation, we only focus on the marginal likelihood functions and make sure that they are accurately approximated.
The particular likelihood function (derived from the modified scenario of coin tossing) properly captures how one parameter interacts with the remaining parameters (as a whole) in the likelihood function.
This is more effective than modeling how one parameter interacts with the remaining parameters individually.
In addition to this particular function, the availability of MLE for parameters as well as the validness of mode mapping makes sure that the marginal likelihood functions are accurate.
More important, the second step of reconstruction is supported by Sklar's Theorem which states that any multivariate joint distribution can be written in terms of univariate marginal distribution functions and a copula which describes the dependence structure between the variables <cit.>. 
Therefore, the second step of reconstruction preserves the precision of the individual marginal likelihood functions, which guarantees that the reconstructed joint likelihood is marginal-wise accurate. 
Third, the two-step approximation is also computationally desirable. 
The major computation cost is from approximating the individual marginal likelihood functions, which is linear in the dimension of the parameter space, making it particularly attractive for large-scale MRFs. 

The rest of the paper is organized as follows. 
Section <ref> systematically presents how we derive this particular likelihood function for marginal likelihood approximation, and how we reconstruct the joint likelihood function via copula.
Section <ref> empirically demonstrates that our MLE-induced likelihood works both more accurate and faster than pseudolikelihood <cit.> and Laplace approximation <cit.>.
Especially as the size of MRFs increases, both the numerical performance and the computational cost of our approach remain consistently satisfactory, whereas the approximation from Laplace approximation deteriorates and pseudolikelihood becomes computationally unbearable. 
We finally conclude in Section <ref>.



§ MLE-INDUCED LIKELIHOOD


In the proposed MLE-induced likelihood, we first approximate the marginal likelihood functions of the individual parameters, and then reconstruct the joint likelihood function using a copula function <cit.>. 
Subsection <ref> introduces a particular parametric likelihood function derived from a modified scenario of coin tossing. 
Subsection <ref> explains how we approximate the marginal likelihood functions with the particular likelihood function, given the availability of the MLE of the parameters. 
Subsection <ref> completes our approximation by linking up the individual marginal likelihood functions via an appropriate copula.



 §.§ A Modified Scenario of Coin Tossing


Suppose that we have a coin with a probability θ (0<θ<1) of landing “head” in a random toss. Let X^i ∈{0,1}, 1 ≤ i ≤ n, denote the landing results, “tails” or “heads”, from n independent tosses of the coin. 
X^i follows a Bernoulli distribution, and the total number of “heads" ∑_i=1^n X^i follows a Binomial distribution. 
We slightly change the conventional setting by casting an external effect, say, a magnetic field, to the coin tossing, so that the probability of landing “head" in a random toss becomes λ, as specified by


    λ = η_0·θ/η_0·θ+(1-η_0)(1-θ),


where η_0 (0<η_0<1) measures the external effect. 
η_0 is multiplicative to θ in the sense given in (<ref>). 
When η_0 is 0.5, this modified scenario reduces to the conventional coin tossing problem.

Under the external effect η_0, suppose that we observe the outcomes from the n random tosses, denoted by 𝕏={x^1,...,x^n}. 
The likelihood function of λ is


    L(λ|𝕏) = P(𝕏 | λ) = nα_0 λ^α_0  (1-λ)^n-α_0,


where α_0 = ∑_i=1^n x^i. Plugging (<ref>) into (<ref>), we have the likelihood function of θ


    L(θ|𝕏) = nα_0    η_0^α_0(1-η_0)^n-α_0
        θ^α_0(1-θ)^n-α_0/(η_0 ·θ + (1-η_0)(1-θ))^n  .


This particular likelihood function (<ref>) is very useful in approximating the marginal likelihood functions of an MRF, as explained next. 



 §.§ Approximating Marginal Likelihood of an MRF


This subsection shows that it is natural to approximate the marginal likelihood functions of an MRF with the particular likelihood function given in (<ref>). 
For illustration, we use a binary pairwise MRF ℳ with d nodes and p parameters. 
The pairwise potential function ϕ_j on edge j (connecting X_u and X_v) parameterized by θ_j (0<θ_j<1) is ϕ_j( X;θ_j)=θ_j^I_{X_u=X_v}(1-θ_j)^I_{X_u ≠X_v}. 
The exact structure of ℳ is provided in Figure <ref>(a).
ℳ includes six nodes and seven edges. 



As shown in Figure <ref>(a), θ parameterizes the edge that directly connects X_1 and X_2. 
The corresponding pairwise potential function is ϕ( X;θ)=θ^I_{X_1=X_2}(1-θ)^I_{X_1≠X_2}. 
Besides this edge, X_1 and X_2 influence each other through the remainder of the network. 
Even for such a small MRF, it is difficult to derive the impact from the remaining part of ℳ towards the concurrence of X_1 and X_2. 
To circumvent the difficulty of explicitly quantifying the overall effect to {X_1,X_2}, we simplify the impact from the remainder of the network regarding the X_1, X_2 concurrence as follows. 

As illustrated in Figure <ref>(b), we mask the remaining part outside {X_1,X_2}.
We assume that the effect from the remaining network can be captured by another potential function. 
Then ℳ is simplified as two nodes X_1 and X_2 and two edges connecting them, one direct edge and one indirect edge masked in the cloud. 
Similar to the parameterization of the original direct edge with θ, we introduce another parameter η to parameterize the potential function on this imaginary indirect edge, which corresponds to the remainder of ℳ. 
From the probability density function of the MRF (a.k.a. a log linear model), θ and η are multiplicative to each other, and the combined effect can be parameterized by  


    λ = η·θ/η·θ+(1-η)(1-θ).


Given the n observed data points 𝕏 from ℳ, the marginal likelihood function of θ can be derived by marginalizing out all of the other parameters. 
By simplifying the remainder of the network into the indirect edge, we lump all of the other parameters into one single parameter η. 




Because it is not feasible in general to obtain the accurate distribution of η, we propose using its point estimate η_0 to approximately quantify its distribution. 

Therefore, we can regard the effect from the remaining network η_0 as the external effect in the modified scenario of coin tossing.
Following the derivation in Subsection <ref>, the marginal likelihood function of θ has a form of (<ref>). 
If the point estimate η_0 is provided, the approximated marginal likelihood function can be easily calculated.

The remaining issue is to find the point estimate η_0.
We rely on the availability of MLE and the fact that “MLE is invariant to reparameterization" <cit.>. 
If the MLEs for θ and λ (denoted by θ̂_MLE and λ̂_MLE) are available, we are able to calculate η_0 by rearranging (<ref>) as


    η_0 = λ̂_MLE· (1-θ̂_MLE)/λ̂_MLE + θ̂_MLE - 2  λ̂_MLE·θ̂_MLE.


In terms of calculation, it is straightforward to show that


    λ̂_MLE = ∑_i=1^n I_{x_1^i=x_2^i}n.


We can get θ̂_MLE from maximum likelihood estimation for an MRF. 
Although the exact likelihood function is not available, a few algorithms <cit.> make use of the concavity of the MRF's log likelihood function and find the MLE via gradient ascent with both satisfactory empirical performance and convergence properties <cit.>.
To find MLE of θ, we use the contrastive divergence algorithm <cit.> or its variations <cit.> which are essentially gradient ascent algorithms. 
Note that the gradient ascent algorithms estimate the MLE from the joint likelihood function. 
Theoretically, the components of the MLE from the joint likelihood function may not necessarily be the same as the MLEs from the marginal likelihood functions, but the difference is expected to be reasonably  small, especially when the joint likelihood function is not highly skewed or multi-modal.

With λ̂_MLE and θ̂_MLE available, η_0 could be calculated by Formula (<ref>).
The marginal likelihood function of θ can be approximated by (<ref>) with α_0 = ∑_i=1^n I_{x_1^i=x_2^i}. 
As one may tell, all the calculation does not depend on the specific structure of the MRF. 
However, our approximation is valid because how one parameter interacts with the remaining parameters in the likelihood function is properly captured.  

Computationally, the approximation of the marginal likelihood functions has no further computational cost other than performing contrastive divergence to find θ̂_MLE. 

Note that this is how we handle θ in the pairwise potential functions on the edges. 
In the situation where θ parameterizes a node potential, the methodology is also applicable. 
In the situation where there are more than two possible outcomes, we need to generalize both Dirichlet distribution and multinomial distribution in a similar fashion.



 §.§ Joint Likelihood Reconstruction From Marginal Likelihood


The marginal likelihood functions provide building blocks for the reconstruction of the joint likelihood function. 
Sklar's Theorem guarantees that any multivariate joint distribution be written in terms of univariate marginal distribution functions and a copula which describes the dependence structure between the variables <cit.>. 
With the marginal likelihood functions available, what we need is a proper copula which specifies the dependence structure between individual parameters <cit.>.

It is rather challenging to capture the accurate copula function given the MRF structure, and the difficulty is expected to grow fast as the dimension of parameter space increases. 
While we prefer a copula function that is as close to the ground truth as possible, it could be impractical to pursue such an overly complicated approximation. 
We aim at a balance between accuracy and feasibility, and consider the conventional copula functions with fewer parameters.

Two types of copula functions are investigated in this work. 
One is the independent copula, and the other one is the exchangeable Gaussian copula. 
The independent copula assumes that the individual marginal likelihood functions in terms of the parameters are independent. 
The joint likelihood function is simply the product of these marginal likelihood functions. 
Namely, if L(θ_j|𝕏) is the marginal likelihood function for θ_j (1 ≤ j ≤ p), then the joint likelihood function is given as follows


    L(θ|𝕏) = ∏_j=1^p L(θ_j|𝕏).


With an exchangeable Gaussian copula, a correlation matrix R needs to be pre-specified. 
R is a p × p positive definite matrix, whose off-diagonal entries are all equal to a certain value ρ <cit.>. 
The copula density in terms of transformed variables is shown below,


    c_R^Gauss( u(θ)) = 1√( R)exp{-12 u(θ)^T (R^-1-I) u(θ)},

where u(θ) = c(u_1,…,u_p) and u_j = Φ^-1 (F_j(θ_j|𝕏)) for 1 ≤ j ≤ p.
Here Φ^-1() is the inverse of cumulative distribution function for the standard normal distribution, and F_j() is the cumulative distribution function from standardizing the marginal likelihood function of θ_j, such as the particular likelihood from the modified scenario of coin tossing. 
When ρ=0, R becomes an identity matrix and this Gaussian copula reduces to the independent copula. 
Given the exchangeable Gaussian copula as well as the marginal likelihood functions, the joint likelihood is specified as

    L(θ|𝕏) = c_R^Gauss( u(θ)) ∏_j=1^p L(θ_j|𝕏).


The simulation results in the next section suggest that the empirical performance of our MLE-induced likelihood is insensitive to the choice of ρ in the exchangeable copula function.



§ SIMULATIONS


In this section, we empirically evaluate our MLE-induced likelihood with two different types of approximation — pseudolikelihood <cit.> and Laplace approximation <cit.>.
Since the likelihood function is heavily used in Bayesian analysis, we evaluate them in MRF Bayesian learning tasks.
Specifically, we compare them in a widely used Markov chain Monte Carlo algorithm — the standard Metropolis-Hastings (MH) algorithm <cit.>.
Note that new types of sampling methods can be used for MRF Bayesian learning, such as Hamiltonian <cit.>, Langevin <cit.>, herding <cit.> and tempered transitions <cit.>, which usually lead to higher efficiency and a better convergence rate.
However since the scope of this paper is to compare different types of approximation of the MRF likelihood function, we do not include these advanced sampling algorithms other than the standard Metropolis-Hastings algorithm.
In Subsection <ref>, we review the setup of the MRF Bayesian parameter learning task and the Metropolis-Hastings algorithm.
In Subsection <ref>, we provide more details about the models and the simulations.
In Subsection <ref>, the empirical results from these different methods are provided and discussed.



 §.§ Bayesian Parameter Learning of MRFs
 

Suppose that an MRF on X∈𝒳^d (𝒳 is a discrete space) is parameterized by θ, and its probability mass function is P( X;θ) =P̃( X;θ) / Z(θ), where P̃( X;θ) is some unnormalized probability measure, and Z(θ) is the normalizing constant, a.k.a. the partition function.
Given a prior of θ and n i.i.d. observed data points 𝕏={ x^1,..., x^n}, Bayesian parameter estimation provides the posterior distribution of θ. This posterior distribution is essential in Bayesian inference. The expected value E(θ|𝕏) (a.k.a. Bayesian estimate or posterior mean) is the optimal estimate under many settings, and the standard deviation of the poster distribution further quantifies the variability of such optimal estimate.

Bayesian parameter estimation for general MRFs is known as doubly-intractable <cit.>. With a prior π(θ) and one data point x, the posterior is P(θ| x) ∝π(θ) P̃( x;θ)/Z(θ). If we use the Metropolis-Hastings (MH) algorithm to generate posterior samples of θ, then in each MH step we have to calculate the MH ratio for the proposal from θ to θ^*


    a(θ^*|θ)     = π(θ^*) P( x;θ^*)Q(θ|θ^*)π(θ) P( x;θ)Q(θ^*|θ)
       =π(θ^*) P̃( x;θ^*)Q(θ|θ^*)Z(θ)π(θ) P̃( x;θ)Q(θ^*|θ)Z(θ^*),


where Q(θ^*|θ) is some proposal distribution from θ to θ^*, and with probability min{1,a(θ^*|θ)} we accept the move from θ to θ^*.
The MH algorithm converges to the correct posterior distribution if the Markov chain is ergodic and in detailed balance.
However, even the MH ratio itself is generally intractable because of the intractable likelihood function.

Primarily, we approximate P( x;θ^*) and P( x;θ) in the MH ratio with our MLE-induced likelihood, pseudolikelihood <cit.> and Laplace approximation <cit.> and compare the performance of the Bayesian estimator when coupled with these three different types of approximation.
We are also aware of other sampling-based methods which make the calculation of the MH ratio feasible via estimating Z(θ)/Z(θ^*), including importance sampling <cit.>, auxiliary variables <cit.>, exchange algorithm <cit.>, and persistent Markov chains <cit.>. 
Although estimating the MH ratio is a broader goal, we still include importance sampling <cit.>, auxiliary variables <cit.>, exchange algorithm <cit.>, and persistent Markov chains <cit.> as secondary baselines in the simulations. 



 §.§ Simulation Setup


We investigate the empirical performance of these different Bayesian estimation methods with three models of different sizes: (i) a 4 × 4 grid MRF with heterogeneous parameters, (ii) a 6 × 6 grid MRF with heterogeneous parameters and (iii) an 8 × 8 grid MRF with heterogeneous parameters.

In total, we compare ten different Bayesian estimation methods, including six likelihood-function-based methods — the exact likelihood calculation (Exact-L), our MLE-induced likelihood from an independent copula (MLE-L, ρ=0), our MLE-induced likelihood from a Gaussian copula with ρ=0.05 (MLE-L, ρ=0.05), our MLE-induced likelihood from a Gaussian copula with ρ=0.1 (MLE-L, ρ=0.1), pseudolikelihood approximation (Pseudo-L), Laplace approximation (Laplace-L), and four sampling-based baseline methods — importance sampling with geometric α (IS-Geometric), auxiliary variable (AuxVar), exchange (Exch) and persistent Markov chain (PersistMC).
For the exact likelihood calculation, we do not enumerate and sum over all possible instantiations for Z(θ), which has a complexity of O(|𝒳|^d).
Instead, we use the efficient recursive algorithm <cit.> which also provides exact calculation of Z(θ) and reduces the computation complexity to O(|𝒳|^l+1), l=min{# rows, # columns} for grid-structured MRFs.
The reason we choose the geometric function for α in the standard importance sampling algorithm is that it is parameter-free and yields accurate estimate of partition function ratios <cit.>.
More algorithmic details of the four secondary baselines are provided in the supplementary material.

In the simulations, we first set the parameters, and then generate a number of data points under the parameters.
We then run the MH algorithm coupled with these different methods of calculating the MH ratio.
When we couple these different methods with the MH algorithm, we set the Markov chains in the same way for a fair comparison, including the initialization of the Markov chain, the number of MH steps and the proposal distribution Q.
For the four sampling-based methods (IS-Geometric, AuxVar, Exch and PersistMC), we use the same number of particles.



 §.§ Simulation Results


In the first set of simulations, we use a 4 × 4 grid-structured MRF with a set of heterogeneous parameters, i.e. θ={θ_1,...,θ_p}, p=24.
The MRF is binary and pairwise. The pairwise potential function on edge j is
θ_j^I_{X_u=X_v}(1-θ_j)^I_{X_u ≠X_v} where nodes X_u and X_v are connected by edge j in the graph.
We set θ_j uniformly distributed on the interval (0.5, 0.8) for j=1,…,p.
We then generate n data points (n=100, 200,..., 1,000).
Eventually, we apply the MH algorithm with ten different methods of calculating the MH ratio, namely Exact-L, MLE-L(ρ=0), MLE-L(ρ=0.05), MLE-L(ρ=0.1), Pseudo-L, Laplace-L, IS-Geometric, AuxVar, Exch and PersistMC. 
We use a uniform distribution for the prior π(θ).
The Q(θ^*|θ) in (<ref>) is a p-variate Gaussian density function with mean θ and covariance matrix σ_Q^2 I_p. The σ_Q^2 is set to be 0.001 in the simulations. 
The total number of MH steps for all methods is set to be 1,000,000.
For the four sampling-based methods (IS-Geometric, AuxVar, Exch and PersistMC), we set the number of particles to be 1,000.
When we have to generate particles from P( X;θ) for a given θ, we randomly instantiate the particles and advance the particles for 1,000 steps under θ.
We observe that coalescence <cit.> always happens within the first 1,000 steps, so the particles can be regarded as exact samples. As one may tell, we set the numbers in a conservative manner, minimizing the impact caused by simulation settings and placing our focus on the genuine performances of methods in comparison.

We use the posterior samples generated from exact likelihood as the benchmark, and compare the posterior samples generated from the other methods by two measures. 
The first is the posterior mean estimated from these methods.
The second is the posterior standard deviation estimated from these methods. 
We keep records of the run time for the methods in the simulations. 
All of the simulations have been replicated 50 times and the averaged results are provided.







Run Time: The run times of these methods in the first set of simulations are listed in Table <ref>.
All the experiments are carried out on one 3GHz Intel Xeon CPU.
Three baselines, IS-Geometric, AuxVar and Exch, did not finish due to their long run time.
We estimate their run time based on how many MH steps they actually finished within the given time.
It is expected that they can finish the 1,000,000 MH steps in thousands of hours.
Persistent Markov chain is the only sampling-based method that finished in a reasonable time (around an hour).
The run time of pseudolikelihood is roughly on the same level, but it is more sensitive to the number of observed data points.
The exact likelihood calculation is tolerable because we use the efficient recursive algorithm <cit.>.
Laplace approximation and our MLE-induced likelihood take the least time, and the run time is insensitive to the sample size.
For the different choices of the copula function, we observe a number of folds increase in run time if a Gaussian copula with a nonzero ρ is selected over the independent copula.

Empirical Performance: Three baselines (IS-Geometric, AuxVar and Exch) did not finish due to their long run time, and therefore their empirical performance cannot be evaluated.
The empirical performance of pseudolikelihood, Laplace approximation, persistent Markov chain and our MLE-L(ρ=0) are plotted in Figure <ref>.
The performance of our MLE-L is insensitive to the choice of the copula function.
Therefore, we do not plot MLE-L(ρ=0.05) and MLE-L(ρ=0.1) in the figure.
For the posterior mean estimated from the posterior samples, it is observed that all of the four methods yield satisfactory results.
As the number of data points increases, the posterior mean estimate from these Bayesian estimators (pseudolikelihood, Laplace approximation, persistent Markov chain and our MLE-L) gets closer to that from exact likelihood estimation.
Persistent Markov chain, Laplace approximation and our MLE-L yield comparable results on this small 4 × 4 grid MRF example.
Pseudolikelihood works slightly worse when there are less data points [The reason is that pseudolikelihood estimators are biased with finite samples, but are asymptotically consistent <cit.>.].
For the estimation of the posterior standard deviation, all the methods provide similar estimates as exact likelihood and the difference decreases as the number of data points increases.
Our MLE-L and persistent Markov chain perform similarly, and yield smaller estimation errors for posterior standard deviation.
Laplace approximation yields a larger estimation error for posterior standard deviation, and pseudolikelihood yields the largest.

Simulations In Higher-dimensional Parameter Space: We replicate the experiments on the 6 × 6 grid-structured MRF with a set of heterogeneous parameters (p=60) and the 8 × 8 grid-structured MRF with a set of heterogeneous parameters (p=112).
We exclude IS-Geometric, AuxVar and Exch due to their long run time.
Exact-L did not finish either, and we lose the benchmark to compare the empirical performance.
Therefore, we compare the Bayesian estimates from pseudolikelihood, Laplace approximation, persistent Markov chain, our MLE-L(ρ=0), MLE-L(ρ=0.05), and MLE-L(ρ=0.1) with the ground truth of the parameters.
All the other configurations in the simulations remain the same.

The run times of the six methods (pseudolikelihood, Laplace approximation, persistent Markov chain, our MLE-L(ρ=0), MLE-L(ρ=0.05), and MLE-L(ρ=0.1)) on the 6 × 6 grid MRF and the 8 × 8 grid MRF are listed in Table <ref>.
Similar to previous simulations on the 4 × 4 grid MRF, Laplace-L and our MLE-L run the fastest, and the run time is insensitive to the number of data points.
For the different choices of the copula function in our MLR-L, there is a number of folds increase in run time if a Gaussian copula with a nonzero ρ is chosen over the independent copula.
The run times of both Pseudo-L and PersistMC are sensitive to the number of data points, and get almost unbearable as both the number of data points increases and the number of parameters increases.

The performance of the six methods (pseudolikelihood, Laplace approximation, persistent Markov chain, our MLE-L(ρ=0), MLE-L(ρ=0.05), and MLE-L(ρ=0.1)) is evaluated by the difference of yielded Bayesian estimates and the ground truth of the parameters.
The numerical results are plotted in Figure <ref>.
The performance of pseudolikelihood, persistent Markov chain and our MLE-L(ρ=0) on the larger MRFs is similar to that from the 4 × 4 grid MRF.
The performance of our MLE-L is also insensitive to the choice of the copula function.
Therefore, we do not plot MLE-L(ρ=0.05) and MLE-L(ρ=0.1) in the figure.
However, the performance of Laplace approximation quickly deteriorates as the dimension of parameter space grows, which suggests Laplace approximation is less favorable than others in large-scale MRFs.







§ DISCUSSION AND CONCLUSION


The primary contribution of our paper is proposing a two-step method of approximating the MRF likelihood function.
The proposed MLE-induced likelihood yields satisfactory numerical precision during the approximation.
Especially as the size of the MRF increases, both the numerical performance and the computational cost of our approach remain consistently satisfactory, whereas Laplace approximation deteriorates and pseudolikelihood becomes computationally unbearable.
We also demonstrate that when our MLE-induced likelihood is used for MRF Bayesian parameter learning, its numerical performance is almost as good as the state-of-the-art sampling-based method — persistent Markov chain <cit.>, but the computational cost is significantly reduced.

From the methodology aspect, our MLE-induced likelihood is desirable for several reasons. 
We break down the approximation into two steps. 
In the first step, we only focus on the marginal likelihood functions and make sure that they are accurately approximated.
The particular likelihood function (derived from the modified scenario of coin tossing) properly captures how one parameter interacts with the remaining parameters (as a whole) in the likelihood function.
This is more effective than modeling how one parameter interacts with the remaining parameters individually (as in Laplace approximation).
During the calculation of the marginal likelihood functions, we also make use of the MLE of the parameters (e.g. from the contrastive divergence algorithm <cit.>), which makes sure that the modes of the marginal likelihood functions are accurate.
The second step of reconstruction is supported by Sklar's Theorem, which preserves the precision of the individual marginal likelihood functions and guarantees that the reconstructed joint likelihood is marginal-wise accurate. 

In addition to the satisfactory empirical performance and the methodology innovation, the proposed MLE-induced likelihood also has some connection with some known distributions.
Especially, the marginal likelihood function we derive in Formula (<ref>) is related to the density function of a Generalized Beta Distribution (GBD) <cit.>.
If the power in the denominator of (<ref>) is n+2 instead of n, function (<ref>) can be standardized to be a density function of a GBD by dividing a constant.
In general, n is relatively large, and n+2 and n are quite close, and (<ref>) is proportional to a GBD density function approximately.

Last but not least, our proposed MLE-induced likelihood nicely extends the recent advances within the field of Markov random field learning. 
Although the typical MLE procedures require that we first specify the log likelihood function, and then find the MLE of the parameters, it does not mean that we cannot find the MLE for MRFs. 
Indeed, a few recent algorithms <cit.> make use of the concavity of the MRF's log likelihood function and find the MLE via gradient ascent with both satisfactory empirical performance <cit.> and convergence properties <cit.>. 
With the availability of MLE of the parameters, we further demonstrate that we are able to recover the intractable likelihood function to some precision from the MLE of the parameters.



plainnat




§ SUPPLEMENTARY: ALGORITHMIC DETAILS OF THE FOUR SECONDARY BASELINES


The real hurdle in Bayesian parameter estimation for general MRFs is the intractable MH ratio 


    a(θ^*|θ)     = π(θ^*) P( x;θ^*)Q(θ|θ^*)π(θ) P( x;θ)Q(θ^*|θ)
       =π(θ^*) P̃( x;θ^*)Q(θ|θ^*)Z(θ)π(θ) P̃( x;θ)Q(θ^*|θ)Z(θ^*).


We can approximate P( x;θ^*) and P( x;θ) in the MH ratio with our MLE-induced likelihood, pseudolikelihood <cit.> and Laplace approximation <cit.>.
Other than approximating the likelihood functions, we can also use other sampling-based methods which make the calculation of the MH ratio feasible via estimating Z(θ)/Z(θ^*), including importance sampling <cit.>, auxiliary variables <cit.>, exchange algorithm <cit.> and persistent Markov chains <cit.>.
Here we provide the algorithmic details of the four secondary baselines.



 §.§ Importance Sampling


We can use importance sampling to estimate r=Z(θ)/Z(θ^*) <cit.> by




    r̂_IS = s_2^-1∑_t=1^s_2P̃( x_2^(t);θ) α( x_2^(t)) s_1^-1∑_t=1^s_1P̃( x_1^(t);θ^*) α( x_1^(t))  ,




where x_1^(1), ..., x_1^(s_1) are sampled from P( X;θ) and x_2^(1), ..., x_2^(s_2) are sampled from P( X;θ^*), and α( X) is an arbitrary function defined on the same support as P̃. 
Theoretically, r̂_IS is a consistent estimator of Z(θ)/Z(θ^*) as long as the sample averages in (<ref>) converge to their corresponding population averages, which is satisfied by Markov chain Monte Carlo under regular conditions. 
However, the optimal choice of α depends on the ground truth of r, and <cit.> provides several options for α, such as a geometric function α( X)=(P̃( X;θ)P̃( X;θ^*))^-1/2.



 §.§ Auxiliary Variables


The second method is to introduce auxiliary variables and cancel Z(θ)/Z(θ^*) in (<ref>). 
<cit.> introduces one auxiliary variable Y on the same space as X, and the state variable is extended to (θ, Y). 
They set the new proposal distribution for the extended state Q(θ, Y|θ^*, Y^*) = Q(θ|θ^*)P̃( Y;θ)/Z(θ) to cancel Z(θ)/Z(θ^*) in (<ref>). 
Therefore by ignoring Y, we can generate the posterior samples of θ via Metropolis-Hastings. 
Technically, this auxiliary variable approach requires perfect sampling <cit.>, but <cit.> pointed out that other simpler Markov chain methods also work with the proviso that they converge adequately to the equilibrium distribution. 
<cit.> extended the single auxiliary variable method to multiple auxiliary variables for improved efficiency, as well as pointed out that the single auxiliary variable method can be simplified as a single-variable exchange algorithm. 
Both the single auxiliary variable algorithm and the single-variable exchange algorithm can be interpreted as importance sampling <cit.>. 
In the auxiliary variable algorithm, r=Z(θ)/Z(θ^*) is estimated by




    r̂_aux = s_1^-1∑_t=1^s_1P̃( x_1^(t);θ̂)P̃( x_1^(t);θ)s_2^-1∑_t=1^s_2P̃( x_2^(t);θ̂)P̃( x_2^(t);θ^*),




where x_1^(1), ..., x_1^(s_1) are sampled from P( X;θ) and x_2^(1), ..., x_2^(s_2) are sampled from P( X;θ^*), and θ̂ is some estimate of θ. 



 §.§ Exchange Algorithm


In the single-variable exchange algorithm, r=Z(θ)/Z(θ^*) is estimated by




    r̂_exch = s^-1∑_t=1^sP̃( x^(t);θ)P̃( x^(t);θ^*),




where x^(1), ..., x^(s) are sampled from P( X;θ^*). 



 §.§ Persistent Markov Chains


Importance sampling, the auxiliary variable method and the exchange algorithm are computationally intensive and do not scale well to large models or high dimensional parameter space. It is because that in each MH step they require generating samples from P( X;θ) for a given θ via the computationally expensive perfect sampling <cit.> or standard Gibbs sampling with long runs. 

In the standard single-variable exchange algorithm, s samples need to be generated from P( X;θ^*) in each MH step when we propose θ^*, and the MH ratio is calculated as (<ref>). The motivation of the persistent Markov chains algorithm <cit.> is that although the proposed θ^* is different from θ, θ^* is usually not far away from θ because MH algorithms usually require proposing small changes so as to maintain a high acceptance rate. Therefore the particles from θ^* should be quite similar to the particles from θ. Therefore, we can reuse the particles generated from θ in the previous step and further advance the particles for k steps under the new parameter θ^*, where k is a small number. Therefore, r=Z(θ)/Z(θ^*) in the l-th MH step can be estimated by




    r̂_persistMC = s^-1∑_t=1^sP̃( x^(t);θ)P̃( x^(t);θ^*),




where x^(1), ..., x^(s) are the particles from the (l-1)-th MH step and get advanced for k steps under θ^*. In our experiments, we set k to be 1. This technique was firstly used in the persistent contrastive divergence algorithm <cit.> which modified the standard contrastive divergence algorithm by reusing the particles persistently to save computation. 


