An Efficient Pseudo-likelihood Method for Sparse Binary Pairwise Markov Network Estimation
Sinong Geng Zhaobin Kuang David Page































||
‖‖
ifstar*
ifstar*





 















prettyitem[1]


































 An Efficient Pseudo-likelihood Method for 
Sparse Binary Pairwise Markov Network Estimation
    Sinong Geng, Zhaobin Kuang, and David Page
 
University of Wisconsin

sgeng2@wisc.edu, zkuang@wisc.edu, page@biostat.wisc.edu
    
=================================================================================================================================




The pseudo-likelihood method <cit.> is one of the most popular algorithms for learning sparse binary pairwise Markov networks. In this paper, we formulate the L_1 regularized pseudo-likelihood problem as a sparse multiple logistic regression problem. In this way, many insights and optimization procedures for sparse logistic regression can be applied to the learning of discrete Markov networks. Specifically, we use the coordinate descent algorithm for generalized linear models with convex penalties <cit.>, combined with strong screening rules <cit.>, to solve the pseudo-likelihood problem with L_1 regularization. Therefore a substantial speedup without losing any accuracy can be achieved. Furthermore, this method is more stable than the node-wise logistic regression approach on unbalanced high-dimensional data when penalized by small regularization parameters. Thorough numerical experiments on simulated data and real world data demonstrate the advantages of the 	proposed method. 




§ INTRODUCTION

Markov networks are a class of probabilistic graphical models that is widely applicable to many areas like image processing <cit.>, multiple testing <cit.> and computational biology <cit.>. In a Markov network, the conditional independence relationships among random variables are illuminated by the structure of the network. The most difficult challenge in estimating binary pairwise Markov networks (BPMNs) has always been dealing with the intractable computation related to log-likelihoods, which makes the learning process an NP-hard problem. Therefore, in the literature, various methods have been proposed to approximate the log-likelihood function instead of using exact estimation. <cit.> built a contrastive divergence algorithm by directly estimating the derivative of log-likelihood for a discrete Markov network. <cit.> considered the neighborhood recovery for each variable separately and proposed the node-wise logistic regression (NLR) method. In <cit.>, pseudo-likelihood (PL) was proposed as an approximation to the log-likelihood of a penalized BPMN. 

According to <cit.>, the PL method is one of the most competitive methods, with faster speed and higher accuracy than the other approaches. However, because of the development in the optimization for other competing methods, the PL method solved by an existing implementation, the  package <cit.>, became the slowest among many learning methods for discrete Markov networks <cit.>; there still remains demand for a more efficient optimization procedure for PL.

Meanwhile, L_1-regularized logistic regression (LR), as one of the most widely used generalized linear models (GLMs), has sophisticated implementations that can deliver solutions efficiently. The state-of-the-art implementation of L_1-regularized GLMs leverages coordinate descent <cit.> with variable screening <cit.>, and has become a building block for many other sparse learning problems like the ones in <cit.>, <cit.> and so on.

As noticed by some researchers, there have been some potential similarities between the objective functions of PL and LR. <cit.> pointed out that the PL model is related to LR problems. <cit.> considered PL as an LR problem with symmetric constraints. <cit.> pointed out the equivalence of the objective functions of the two problems and provided an optimization algorithm based on this relationship. However the advantages of treating the sparse PL model as a regularized LR problem has been neither fully exploited nor emphasized enough. As a result, in a recent empirical study that compares multiple estimation mehtods for BPMNs <cit.>, optimization procedures with suboptimal efficiency are still considered and benchmarked to solve PL problems.


By posing PL as an LR problem in the context of learning an L_1-regularized BPMN, a much faster alternative to the optimization of PL can be attained. Specifically our work in this paper is summarized as follows:




  * With the relationship between PL and LR, we consider an optimization procedure using the coordinate descent algorithm and variable screening to solve the L_1 regularized PL problem.  The procedure in question can be conveniently implemented via the state-of-the-art optimization algorithm for learning sparse generalized linear models: . Thus, the procedure is called PLG (pseudo-likelihood using glmnet). Achieving a dramatic speedup without losing any accuracy, PLG substantially outperforms the highly visible implementation of PL <cit.>. 

  * Unlike the NLR approach, the PLG procedure maintains the efficiency even when dealing with unbalanced high-dimensional data penalized by small regularization parameters. We provide insights and numerical experiments to demonstrate the superior stability of the PLG method.     
	   



§ BACKGROUND


To motivate the optimization approach for PL, we first review the background knowledge about BPMNs. 

We consider a p-dimensional binary observation = (x_1, x_2, …, x_p)^⊤∈{0,1}^p. In a BPMN, the distribution of  is associated with a network with vertex set V = {1, 2, …, p} and edge set E ⊆ V × V. Accordingly, based on the ground truth parameter Θ^*, a p × p symmetric matrix given as: 

    Θ^* = 
    [ θ^*_11 θ^*_12      … θ^*_1p; θ^*_21 θ^*_22      … θ^*_2p;      ⋮      ⋮      ⋱      ⋮; θ^*_p1 θ^*_p2      … θ^*_pp ]
     = [ θ^*_ij]_p × p,

the joint probability mass function  is defined as:

    P_Θ^*() = exp(∑_s∈ Vθ^*_ss x_s + ∑_(s,t) ∈ Eθ^*_st x_s x_t - Ψ(Θ^*)),

where Ψ(Θ), for any symmetric Θ = [ θ_ij]_p × p, denotes the log-partition function defined as:

    Ψ(Θ) = log[∑_∈{0,1}^p exp(∑_t⩾ s ⩾ 1θ_st x_s x_t ) ].

 
In order to estimate (<ref>), people consider the L_1-regularized log-likelihood function for the BPMN:

    ℒ(Θ ;) = ∑^p_t⩾ s ⩾ 1θ_st (^⊤)_st - NΨ(Θ) - N λ/2∑_s ≠ tθ_st,
 
where (^⊤)_st denotes the element in the s^th row and the t^th column of , given N independent and identically distributed samples = (_1, _2, …, _N)^⊤ = [ _ij]_N × p. The goal of dealing with (<ref>) is to estimate Θ^* with the optimal Θ given as: 

    Θ_^* = _Θℒ(Θ,),

which can be can be extremely challenging because of the intractable log-partition function, Ψ(Θ). Therefore, instead of maximizing (<ref>), the PL method considers the L_1-regularized pseudo-likelihood function <cit.>,

    ℒ̂(Θ;) = 
    ∑_n=1^N 
    ∑^p_s = 1[_ns( θ_ss + ∑_t≠ s_ntθ_st) 
    
    - Ψ_s(_n, Θ) ]
    - N λ∑_t > sθ_st,

as a replacement for the penalized log-likelihood function (<ref>) and solves for 

    Θ̂_^* = _Θℒ̂(Θ,)

to estimate Θ^*.
Here, Ψ(Θ) is replaced by the much simpler 

    Ψ_s(; Θ) = log[1 + exp(θ_ss + ∑_t ≠ s x_t θ_st)].
  

Compared with exact methods that solve (<ref>), the PL method that solves (<ref>) is shown to be more efficient without sacrificing too much accuracy <cit.>.

There is also a connection between (<ref>) and the objective function of the NLR algorithm, which separately maximizes 

    ∑_n=1^N
    [_ns( θ_ss + ∑_t≠ s_ntθ_st)-  Ψ_s(_n, Θ) ]
    
    - N λ∑_t > sθ_st,

for all s ∈ V. In fact, (<ref>) is just the sum of (<ref>) on s's. Specifically, <cit.> considered maximizing (<ref>) as an L_1-regularized LR problem with response 

    y_s = (_1,s _2,s … _N,s)^⊤.




§ CONVERSION FROM A PSEUDO-LIKELIHOOD PROBLEM TO A SPARSE LOGISTIC REGRESSION PROBLEM


We now demonstrate the relationship between the objective functions in PL and LR by transforming (<ref>) into a logistic loss function with parameter , design matrix , and response , which are defined subsequently.

We first define parameter . Since in LR problems the parameter is a vector instead of a matrix like Θ, we redefine the parameter in BPMNs by stacking the upper triangular elements of Θ column by column to a vector and appending the diagonal elements to the end. Thus we have


    =
     [ upper triangular elementsθ_12,     θ_13,    θ_23, … , θ_(p-1)p,                             diagnoal elementsθ_11,   …,    θ_pp;                                                                 ],
   
where for any (s,t) ∈{(s,t)|s ≠ t,  (s,t) ∈ V × V }, θ_st is transformed into the j_st^th element in  with 

    j_st = min(s,t)+(max(s,t)-2)(max(s,t)-1)/2.
 

Then for the definition of matrix = [ _ij]_Np × (m+p), with m = p(p-1)/2 denoting the number of upper triangular elements, we review the concept of the indicator function: 

    1( C ) = 
    1    C  is sataified
    
    0   otherwise
    .

Furthermore, we define 

    _ij = _n t   ∃  (s,t)  s.t.  j = j_st
    1 (s=j-m)     m+1 ⩽ j ⩽ m+p
    
    0    otherwise,

where i = N(s-1)+n, n ∈{1,2, …, N}. 

Finally, the response  is defined as:

    = (_11, _21, …, _N1, …, _1p, _2p,…, _N p)^⊤.


With , , and , we can rewrite the first part (the log-likelihood) of (<ref>) as:

    ^⊤ - ∑_n = 1^N∑_k=1^pΨ_k(_n, Θ)
    
    =    ^⊤ -  ∑_k=1^Nplog[1 + exp(()_k^⊤)],

where

    ()_k = (_k1, _k2, …, _k(m+p))^⊤,

and

    ∑_n = 1^N∑_k=1^pΨ_k(_n, Θ) = ∑_k=1^Nplog[1 + exp(()_k^⊤)]

because of (<ref>). Therefore, (<ref>) is equal to

    ^⊤ - ∑_k=1^Nplog[1 + exp(()_k^⊤)] -N λ∑_s ≠ tθ_st,
 
which is exactly the loss function for a penalized LR problem consisting of Np samples with the design matrix , response , and the parameter . 

As a result, we have converted an L_1-regularized pseudo-likelihood problem into an LR problem with the objective function (<ref>). 

Based on the relationship established above, we can solve a spare PL problem by solving its equivalent sparse LR problem. The consequence is that we can take advantage of the sophisticated optimization procedures for sparse LR problems to compute the solution for a PL problem efficiently. We now consider the optimization procedure for sparse PL problems based on one of the most efficient optimization algorithms for penalized LR problems, the coordinate descent algorithm <cit.>. Furthermore, we also use the initialization procedure provided by the strong screening rule <cit.> for a further speedup. The details of PLG are presented in Algorithm <ref>. 

[algorithm]font=footnotesize


 




§ PLG VERSUS NLR WHEN LEARNING FROM UNBALANCED DATA


In this section, we compare the performance of the PLG method with the NLR algorithm when dealing with unbalanced high-dimensional data under small regularization parameters.

When the data are unbalanced, some columns of  will be predominantly 0's or 1's. Since NLR uses each column of  in turn as a response, columns with predominantly 0's or 1's will also serve as responses in some spare LR models. Unfortunately, in this situation, we observe that NLR will be extremely slow and even may fail to converge. To make things worse, since the optimal solution gets denser with the decrease of the regularization parameter, the efficiency of NLR further deteriorates when penalized by small λ's.

On the contrary, PLG is more stable in this situation. In the PLG method, the response is chosen to be a longer vector  in (<ref>), consisting of all the elements in  instead of just one column. Therefore, the unbalanced nature of one column in  cannot have a huge effect on . That is to say, while NLR is dealing with highly unbalanced responses, PLG is still solving an ordinary LR problem with a relatively balanced response, and thus maintains high efficiency even in the context of unbalanced high-dimensional data penalized by small regularization parameters. We will illustrate this phenomenon in detail in Section <ref> by experiments. 





§ NUMERICAL EXPERIMENTS

In this section, we compare the empirical performance of the PLG algorithm with those of the original implementation of the PL method by  (BMN) <cit.> and the NLR algorithm <cit.>. First, we show that, as a more efficient optimization method, PLG achieves the same objective function value as BMN with a much faster speed. We also apply NLR and PLG to unbalanced high-dimensional data for a deeper understanding of the relative efficiency of the two methods. Furthermore, we evaluate the structure learning performances  of the three methods using  receiver operating characteristic (ROC) curves as our evaluation metric. ROC curves of the three methods are shown to be very similar on simulated datasets. Finally, performances of the three methods on real world data are also investigated using the senator voting record <cit.> dataset.




Except BMN, the existing best implementation of PL methods, we are only comparing PLG with NLR, while many estimation methods for sparse BPMNs have been proposed in the literature <cit.>. We believe that the comparison between the performance of PLG and that of NLR is representative because NLR has been empirically shown to have the highest efficiency compared with other competing methods especially when dealing with high-dimensional data <cit.>. Furthermore, as shown in Section <ref>, NLR also uses a kind of pseudo-likelihood to approximate log-likelihood and thus has a similar objective function to that of PLG, making the contrast between PLG and NLR a natural comparison.   



 §.§ Simulated Data Generation
  
We use a procedure similar to that in  <cit.> to generate the ground truth parameter Θ^* and the synthetic datasets.
*

  * The number of the vertices, p, in the ground truth network, is chosen to be 5, 10, 15, 20 or 25. Each element of Θ is drawn randomly to be non-zero, with edge generation probability 𝐏∈{0.2, 0.3, 0.4, 0.5}. And the non-zero elements have a uniform distribution on [-1,1].

  * 1000 samples are generated by Gibbs sampling with 1000 burn-in steps.

  * The results reported from Section <ref> to Section <ref> are averages of 20 trials.





 §.§ Model Selection


Before we proceed to compare the efficiency and accuracy of BMN, NLR, and PLG, we conduct model selection to find the best regularization parameters (best representatives) for the three methods in different networks. We use StARS, a stability-based regularization parameter selection method to high dimensional inference for undirected graphs <cit.>, to determine the λ that achieves the best balance between the edge selection stability and the network sparsity. In addition, since BMN and PLG have the same objective function, it's reasonable to use the same λ's for them and the λ's for NLR should be half of those in BMN <cit.>. In detail, for different networks and methods, the λ's selected are summaried in Table <ref>. 






 §.§ Efficiency
  

We now compare the efficiency among BMN, NLR and PLG. Using the λ's chosen by the model selection procedure in Section <ref>, we apply the three methods to datasets generated by BPMNs with different numbers of vertices and edge generation probability 𝐏s. The computation time of the learning process after the selection of λ's is reported in Figure <ref>. 

With the improvement on optimization, PLG outperforms BMN tremendously and becomes comparable to NLR in the aspect of efficiency. Furthermore, the advantages of NLR and PLG to BMN are more substantial with the scaling up of networks. This observation is consistent with the existing result that NLR performs better for high-dimensional problems <cit.>. 

In fact, the comparable performances between NLR and PLG are also reasonable, considering that they apply similar optimization approaches to similar objective functions. 

Naturally, we want to see whether PLG compromises accuracy for acceleration. To this end, we examine the difference between the parameters achieved by PLG (_P) and that achieved by BMN (_B). For a clear illustration of the difference, we define the relative difference ϵ between _P and _B as

    ϵ = _P-_B_2/_B_2.

The ϵ's for the solutions achieved by PLG and BMN in the experiments above are shown in Table (<ref>). It should be noticed that PLG achieves nearly the same parameters as BMN. In addition, the relative difference rises with the increase of p since we are using the same stopping criterion in different networks. These results suggest that PLG provides a notable improvement in efficiency to PL methods without losing any accuracy.   







 §.§ Structure Estimation
  
In order to demonstrate that PLG also has the ability to estimate the correct structure of ground truth networks we compare ROC curves of the three methods for structure estimation to networks consisting of 15 vertices with edge generation probability 𝐏∈{ 0.2, 0.3, 0.4, 0.5} respectively in Figure (<ref>). Overall, we notice that all the three methods achieve nearly the same performance, and in some figures, the lines even overlap with each other, demonstrating the utility of PLG in structure estimation. This result is not surprising because of the coherence with existing empirical results in the literature <cit.>. 




 §.§ Learning from Unbalanced Data

 

As we mention above, with the PLG implementation, PL becomes comparable to NLR in efficiency and accuracy. In fact, PLG even outperforms NLR when dealing with unbalanced data under small λ's. 

Analysis on this particular situation is necessary because using the joint probability mass function (<ref>), it is very likely to generate extremely unbalanced data if θ_ss≠ 0 in a BPMN. In addition, although we only use simulated data to illustrate the stability of two methods for unbalanced data with small λ's, unbalanced data are ubiquitous in practical problems <cit.>, especially in researches on mutations of genes. Considering the huge amount of genes and the small probability of mutations, only a few mutated samples can be observed in practice. Accordingly, samples in this kind of problems are always unbalanced and thus it's meaningful to scrutinize the extreme situation in our work.

To compare the stability of PLG and NLR with unbalanced data and small λ's, we set the number of the vertices, p=10 and compare the computation time under different λ's. In addition, we assume θ_1,1=5, when simulating samples to generate unbalanced data. We omit the results of BMN because of the low efficiency. Results are summarized in Figure <ref>.

We notice that the efficiency of NLR decays immediately with small λ's and unbalanced data while our method still maintains a fast speed, as we expect in section <ref>. PLG performs similarly with NLR under large λ's, which is consistent with the results in <ref>. Theses results indicate that PLG is readily available for unbalanced data. 




 §.§ Real World Experiments

In this section, performances of the three methods on real world data are investigated. 
We conduct an experiment using the senator voting record consisting of 279 samples and 100 variables in the second session of the 109th Congress<cit.>. The task of interest is to investigate the clustering effect of voting consistency. That is to say, we want to find the senators who are more likely to cast similar votes on bills. Here are the details of the experiment. 
*

  * Each bill is considered as a sample and the votes from senators are features.

  * If a senator votes for one bill, the corresponding element in the sample will be denoted by 1, otherwise 0. Missing data are imputed as 0.

  * A binary pairwise Markov network is used to model the data. And we learn the network with PLG and BMN.

  * The λ selected by the model selection procedure in Section <ref> is 0.06.  

  * The vertices represent senators and the edges denote the estimated θ_st,where s≠ t. Furthermore, only the edges with a positive estimated parameter are displayed.

The visualization of the BPMN for senator voting data is presented in Figure <ref>. First, as would be expected, senators are divided into mainly a Democratic cluster and a Republican cluster, which is roughly consistent with party memberships of the senators. Second, as a Republican, Lincoln Chafee is closer and has more connection to Democrats. In fact, he joined the Democrats in 2008 <cit.> and our model can detect his democratic-leaning voting pattern. Third, since Ben Nelson is “one of the most conservative Democrats" <cit.>, it's not surprising that his voting record is closer to Republicans and has some connections with them. 

Furthermore, senators in the same parties and the same states tend to have more connections. These findings all coincide with conventional wisdom, suggesting that PLG can capture interesting dependencies in practical problems. 
We also contrast the efficiency of NLR and PLG on this real world data under different λ's in Figure <ref>. Again, the results of BMN are not included because of its long runtime. Similar to the results in Figure <ref>, the computation time is almost the same for the two methods, indicating the high efficiency of PLG for practical problems.



§ CONCLUSION

For the task of accelerating the optimization of PL, by the equivalence between the objective functions of PL and LR we have studied an optimization method for PL models in the context of BPMNs. Experimental results suggest that PLG is a viable candidate towards scalable and efficient learning of BPMNs even in extreme conditions. Although we focus on binary pairwise Markov networks, our method is generally applicable to other discrete Markov networks whose pseudo-likelihood functions have a close relationship to logistic loss functions.  






named



