High-dimensional structure learning of binary pairwise Markov networks: A comparative numerical study

TL;DR

This paper compares various methods for learning the structure of high-dimensional binary Markov networks, demonstrating that pairwise tests can outperform pseudo-likelihood approaches in accuracy for certain applications.

Contribution

It provides a comprehensive numerical comparison of pseudo-likelihood and pairwise test methods for high-dimensional binary Markov network structure learning.

Findings

Pairwise methods can be more accurate than pseudo-likelihood methods in high-dimensional settings.

A parallel Gibbs sampler based on restricted Boltzmann machines is proposed for efficient sampling.

The study offers insights into method performance across different high-dimensional scenarios.

Abstract

Learning the undirected graph structure of a Markov network from data is a problem that has received a lot of attention during the last few decades. As a result of the general applicability of the model class, a myriad of methods have been developed in parallel in several research fields. Recently, as the size of the considered systems has increased, the focus of new methods has been shifted towards the high-dimensional domain. In particular, introduction of the pseudo-likelihood function has pushed the limits of score-based methods which were originally based on the likelihood function. At the same time, methods based on simple pairwise tests have been developed to meet the challenges arising from increasingly large data sets in computational biology. Apart from being applicable to high-dimensional problems, methods based on the pseudo-likelihood and pairwise tests are fundamentally very different. To compare the accuracy of the different types of methods, an extensive numerical study is performed on data generated by binary pairwise Markov networks. A parallelizable Gibbs sampler, based on restricted Boltzmann machines, is proposed as a tool to efficiently sample from sparse high-dimensional networks. The results of the study show that pairwise methods can be more accurate than pseudo-likelihood methods in settings often encountered in high-dimensional structure learning applications.

Figures

Tables

Table 1. Table 1: Various structural characteristics of the generated network structures ( mean ± sd plus-or-minus mean sd \text{mean}\pm\text{sd} ). The average/maximum neighbourhood size is the average/maximum number of nodes connected to a node. The average clustering coefficient measures the average connectedness between neighbours of a node.

Network property	$d$	Grid network	Scale-free network	Small-world network
	40	$3.35\pm 0$	$4.00\pm 0.09$	$4\pm 0$
Avg. neigbourhood size	200	$3.70\pm 0$	$4.15\pm 0.04$	$4\pm 0$
	1000	$3.91\pm 0$	$4.20\pm 0.02$	$4\pm 0$
	40	$4\pm 0$	$10.16\pm 1.17$	$6.24\pm 0.69$
Max. neigbourhood size	200	$4\pm 0$	$14.40\pm 1.60$	$6.84\pm 0.58$
	1000	$4\pm 0$	$19.40\pm 2.08$	$7.82\pm 0.60$
	40	$0\pm 0$	$0.11\pm 0.03$	$0.27\pm 0.05$
Avg. clustering coefficient	200	$0\pm 0$	$0.03\pm 0.01$	$0.23\pm 0.02$
	1000	$0\pm 0$	$0.01\pm 0.00$	$0.23\pm 0.01$

Table 2. Table 2: Average runtime of methods (in seconds). The table shows the runtime of the postprocessing step of the pairwise methods. The total runtime of a pairwise method is obtained by adding the runtime of MI.

Method	$d$	Runtime
		$n=200$	$n=1000$	$n=5000$
	40	0.06	0.02	0.07
MI	200	0.40	0.53	1.34
	1000	8.99	12.38	28.50
	40	0.05	0.00	0.00
CLR	200	0.03	0.02	0.02
	1000	2.35	2.11	2.01
	40	0.22	0.18	0.18
ARACNE	200	5.12	5.12	5.16
	1000	173.73	176.73	177.69
	40	0.13	0.00	0.00
ND	200	0.06	0.03	0.03
	1000	2.39	2.44	2.53
	40	17.07	18.28	37.18
PlmL1	200	178.55	242.67	594.57
	1000	1729.01	8028.80	31062.33
	40	0.62	1.62	6.42
PlmDCA	200	10.25	37.50	173.20
	1000	240.88	1286.27	8171.49

Equations

$$p(x;\mathbf{\theta })=\frac{1}{Z(\mathbf{\theta })}\exp \left(\sum _{v\in V}{x}_v{\theta }_v+\sum _{\{v,v^{\prime }\}\in E}{x}_v{x}_{v^{\prime }}{\theta }_{v{v}^{\prime }}\right),$$

$$Z(\mathbf{\theta })=\sum _{x\in \mathcal{X}}\exp \left(\sum _{v\in V}{x}_v{\theta }_v+\sum _{\{v,v^{\prime }\}\in E}{x}_v{x}_{v^{\prime }}{\theta }_{v{v}^{\prime }}\right).$$

$$p(x;A,b)=\frac{1}{Z(A,b)}\exp \left(\frac{1}{2}{x}^{\top }Ax+b^{\top }x\right),$$

$$p(x,y;W,b)\propto \exp \left(y^{\top }Wx-\frac{1}{2}{y}^{\top }y+b^{\top }x-\frac{1}{2}\text{diag}(W^{\top }W{)}^{\top }x\right),$$

$$\{\begin{array}{l}{W}_{(v{v}^{\prime },v)}=\sqrt{|{\theta }_{v{v}^{\prime }}|}\\ W_{(v{v}^{\prime },v^{\prime })}=\text{sign}({\theta }_{v{v}^{\prime }})\sqrt{|{\theta }_{v{v}^{\prime }}|}\end{array}\text{ for all }\{v,v^{\prime }\}\in E.$$

$$p(x;W,b)\propto \exp \left(\frac{1}{2}{x}^{\top }{W}^{\top }Wx+b^{\top }x-\frac{1}{2}\text{diag}(W^{\top }W{)}^{\top }x\right).$$

$$\frac{1}{2}{x}^{\top }{W}^{\top }Wx=\frac{1}{2}{x}^{\top }Ax+\frac{1}{2}\text{diag}(W^{\top }W{)}^{\top }x,$$

$$A=\left(\begin{array}{ccc}0& {\theta }_{12}& 0\\ {\theta }_{12}& 0& -{\theta }_{23}\\ 0& -{\theta }_{23}& 0\end{array}\right),$$

$$W=\left(\begin{array}{ccc}\sqrt{{\theta }_{12}}& \sqrt{{\theta }_{12}}& 0\\ 0& \sqrt{{\theta }_{23}}& -\sqrt{{\theta }_{23}}\end{array}\right)\text{ and }{W}^{\top }W=\left(\begin{array}{ccc}{\theta }_{12}& {\theta }_{12}& 0\\ {\theta }_{12}& {\theta }_{12}+{\theta }_{23}& -{\theta }_{23}\\ 0& -{\theta }_{23}& {\theta }_{23}\end{array}\right).$$

$$Y_{v{v}^{\prime }}\mid X=x\sim N(W_{(v{v}^{\prime },:)}x,1)\text{ for }\{v,v^{\prime }\}\in E.$$

$$p(X_v=x_v\mid Y=y)=\sigma \left(({\theta }_v-\frac{1}{2}(W^{\top }W{)}_{(v,v)})x_v+y^{\top }{W}_{(:,v)}\right)\text{ for }v\in V,$$

$$l(\mathbf{\theta };\mathbf{x})=\log p(\mathbf{x};\mathbf{\theta })=\sum _{i=1}^n\log p(x^{(i)};\mathbf{\theta }).$$

$$pl(\mathbf{\theta };\mathbf{x})=\sum _{v=1}^d\log p({\mathbf{x}}_v\mid {\mathbf{x}}_{-v};{\mathbf{\theta }}_v)=\sum _{v=1}^d\sum _{i=1}^n\log p(x_v^{(i)}\mid x_{-v}^{(i)};{\mathbf{\theta }}_v),$$

$$p(x_v\mid x_{-v};{\mathbf{\theta }}_v)=\frac{1}{Z_v({\mathbf{\theta }}_v)}\exp \left(x_v{\theta }_v+\sum _{v^{\prime }\in -v}{x}_v{x}_{v^{\prime }}{\theta }_{v{v}^{\prime }}\right),$$

$$Z_v({\mathbf{\theta }}_v)=1+\exp \left({\theta }_v+\sum _{v^{\prime }\in -v}{x}_{v^{\prime }}{\theta }_{v{v}^{\prime }}\right)$$

$$\underset{\mathbf{\theta }}{\mathrm{argmax}}\sum _{v=1}^d((\mathbf{x}{\mathbf{\theta }}_v{)}^{\top }{\mathbf{x}}_v-n\log Z_v({\mathbf{\theta }}_v)).$$

$$\underset{{\mathbf{\theta }}_v}{\mathrm{argmax}}(\mathbf{x}{\mathbf{\theta }}_v{)}^{\top }{\mathbf{x}}_v-n\log Z_v({\mathbf{\theta }}_v)v=1,\dots ,d.$$

$$\underset{{\mathbf{\theta }}_v}{\mathrm{argmax}}(\mathbf{x}{\mathbf{\theta }}_v{)}^{\top }{\mathbf{x}}_v-n\log Z_v({\mathbf{\theta }}_v)-\lambda \Vert {\mathbf{\theta }}_v{\Vert }_{1}v=1,\dots ,d,$$

$$\min \left(|{\theta }_{v{v}^{\prime }}|,|{\theta }_{v^{\prime }v}|\right)>0,$$

$$\max \left(|{\theta }_{v{v}^{\prime }}|,|{\theta }_{v^{\prime }v}|\right)>0.$$

$$\underset{{\mathbf{\theta }}_v}{\mathrm{argmax}}(\mathbf{x}{\mathbf{\theta }}_v{)}^{\top }{\mathbf{x}}_v-n\log Z_v({\mathbf{\theta }}_v)-\lambda \Vert {\mathbf{\theta }}_v{\Vert }_2^{2}v=1,\dots ,d,$$

$$MI(v,v^{\prime })=\sum _{x_v\in {\mathcal{X}}_v}\sum _{x_{v^{\prime }}\in {\mathcal{X}}_{v^{\prime }}}p(x_v,x_{v^{\prime }})\log \left(\frac{p(x_v,x_{v^{\prime }})}{p(x_v)p(x_{v^{\prime }})}\right).$$

$$\hat{p}(x_v,x_{v^{\prime }})=\frac{n(x_v,x_{v^{\prime }})+1/4}{n+1}\text{ and }\hat{p}(x_v)=\sum _{x_{v^{\prime }}\in {\mathcal{X}}_{v^{\prime }}}\hat{p}(x_v,x_{v^{\prime }}),$$

$$f(v,v^{\prime })=\sqrt{z_v(v,v^{\prime }{)}^2+z_{v^{\prime }}(v,v^{\prime }{)}^2}.$$

$$MI(v,v^{\prime })\le \min \left(MI(v,v^{\prime \prime }),MI(v^{\prime },v^{\prime \prime })\right).$$

$$G_{\text{obs}}=G_{\text{dir}}+G_{\text{indir}}=G_{\text{dir}}+(G_{\text{dir}}^2+G_{\text{dir}}^3+\dots ).$$

$$G_{\text{dir}}=G_{\text{obs}}(I+G_{\text{obs}}{)}^{-1}.$$

$$|{\theta }_v|\sim \text{Uniform}(0,1)\text{and}|{\theta }_{v{v}^{\prime }}|\sim \text{Uniform}(1,2),$$

$$\text{precision}=\frac{TP}{TP+FP}\text{and}\text{recall}=\frac{TP}{TP+FN}.$$

$$\lambda \in \{\frac{k}{100}{\lambda }_{\max }{\}}_{k=1}^{100},$$

Full text

1 Introduction

Learning the dependency pattern over a large collection of variables is an important problem encountered in various fields of science [1]. This learning task can be cast as a sub-problem in the popular class of probabilistic graphical models. The ultimate goal of graphical models is to represent the joint distribution over a collection of variables by exploiting statements of conditional independence. The key feature of graphical models is to use a graph structure to compactly encode the dependence structure over the variables. In addition to being an essential part of graphical models, the graph structure provides a natural target for applications where the main goal is to gain an understanding in how variables within a system interact with each other.

In this work, we will consider the problem of learning the undirected graph structure of pairwise Markov networks over binary variables. As a result of the general applicability of Markov networks, new and exciting structure learning methods have been developed in parallel in various fields including machine learning [2], statistics [3], statistical physics [4, 5], and computational biology [6, 7, 8, 9]. In particular, as the size of the considered problems has increased, computational scalability of the algorithms has become more important. For example, the huge number of variables encountered in genome-wide sequencing studies in computational biology has steered the research towards methods based on simple pairwise tests [6, 7, 8, 9]. Although pairwise methods do not come with any general asymptotic guarantees that would hold in a general setting, they have been shown to perform well on various real-world applications. In contrast, the pseudo-likelihood methods, which have been developed in the more theoretical fields of statistics [3] and statistical physics [4, 5], are more elaborate and enjoy nice theoretical properties, such as consistency.

Although pairwise and pseudo-likelihood methods have been designed for the same underlying problem, they have rarely been compared against each other in a controlled setting. In particular, the finite-sample performance of the pseudo-likelihood methods in comparison to the pairwise methods has not been thoroughly examined. Therefore, in this work, we perform a numerical simulation study on synthetic models to compare a collection of methods from the two algorithmic families. The focus of the study is on high-dimensional sparse models, which are often encountered in real-world applications. To facilitate sampling from such models, we present a parallelizable and highly scalable Gibbs sampler based on sparse restricted Boltzmann machines.

We begin in Section 2 by introducing binary pairwise Markov networks. After that, we continue by presenting the Gibbs sampler and introducing the structure learning problem. In Section 3, we present the methods included in the study. We discuss the main idea and the underlying assumptions behind each method. In Section 4, we go through the experimental setup and present the key findings. We conclude this work with a discussion in Section 5.

2 Pairwise Markov networks

A Markov network (or Markov random field) is an undirected graphical model that represents a joint distribution over a set of $d$ random variables $X=\{X_{1},\ldots,X_{d}\}$. In this work, we assume a positive distribution and we also restrict the variables to be binary with outcome space $\mathcal{X}_{v}=\{0,1\}$. We use a lowercase letter $x_{v}$ to denote that the corresponding variable has been assigned a specific value in $\mathcal{X}_{v}$. The dependence structure of a Markov network is encoded by an undirected graph $G=(V,E)$, where $V=\{1,\ldots,d\}$ are nodes (or vertices), corresponding to the variables, and $E\subseteq V\times V$ are edges, representing direct dependencies between the variables. In particular, absence of edges implies statements of conditional independence, such that separation between nodes in the graph implies conditional independence between the corresponding variables (see e.g. [1] for more details). The conditional independence statements are the fundamental assumptions used by graphical models to decompose the joint distribution into smaller components.

In this work, we will consider a special subclass of Markov network that only includes pairwise interactions. The (positive) joint distribution over the variables can then be parameterized as the log-linear model:

[TABLE]

where $\boldsymbol{\theta}$ contains the model parameters, which consist of the node-specific bias parameters $\{\theta_{v}\}_{v\in V}$ and the edge-specific interaction parameters $\{\theta_{vv^{\prime}}\}_{\{v,v^{\prime}\}\in E}$ [10]. The function $Z(\boldsymbol{\theta})$ is a normalizing constant known as the partition function:

[TABLE]

For notational convenience, both $\theta_{vv^{\prime}}$ and $\theta_{v^{\prime}v}$ denote the same interaction parameter, unless otherwise stated. Binary pairwise Markov networks are equivalent to Ising models, which are a well-studied model class in the statistical physics community.

2.1 Sampling

Sampling from a non-chordal Markov network is an important and non-trivial problem. The common approach is to use MCMC methods, the most well-known being the sequential Gibbs sampler. In each iteration, the algorithm runs through the variables in a pre-specified order and updates the value of each variable by conditioning on the rest of the variables or, for a specified graph, the neighbors of the variable. In this work, we will use an alternative MCMC sampler based on sparse Restricted Boltzmann Machines (sRBM) with real-valued hidden variables. The approach is similar to method in [11], where an RBM is used to sample fully-connected networks.

Let $G=(V,E)$ be the graph of a binary pairwise Markov network. Also, for the sake of the derivation in this section, consider the following matrix-based representation of the distribution:

[TABLE]

where $A$ is an interaction matrix, such that $A_{(v,v^{\prime})}=\theta_{vv^{\prime}}\text{ if }\{v,v^{\prime}\}\in E$ and 0 otherwise, and $b=(\theta_{1},\ldots,\theta_{d})$ is a vector containing the bias parameters.

For each $\{v,v^{\prime}\}\in E$, we introduce an auxiliary continuous variable, which we (with some abuse of notation) denote by $Y_{vv^{\prime}}$. Thus, the number of auxiliary variables will be equal to the number of edges, which is denoted by $n_{\text{edges}}$. Now, the joint distribution over $(X,Y)$ will be defined in form of an sRBM, where the original variables $X$ are visible and the auxiliary variables $Y$ are hidden:

[TABLE]

where $W$ is a sparse $n_{\text{edges}}$-by-$d$ edge weight matrix connecting $X$ and $Y$. Note that there are no direct interactions within $X$ or $Y$. Let $W_{(vv^{\prime},:)}$ represent the row corresponding to variable $Y_{vv^{\prime}}$ and assume that $v<v^{\prime}$. Each row $W_{(vv^{\prime},:)}$ will contain two non-zero elements in columns $v$ and $v^{\prime}$, which connect $X_{v}$ and $X_{v^{\prime}}$ indirectly via $Y_{vv^{\prime}}$. More specifically, to maintain the original edge-specific interactions between the variables in $X$ in the sRBM, the non-zero elements in $W$ are specified according to

[TABLE]

Following the same line of reasoning as [11], the marginal distribution over $X$ can be shown to be

[TABLE]

As a result of how $W$ was constructed, we further have that

[TABLE]

which ultimately means that the marginal distribution (2) is equivalent to (1). In other words, we have constructed an sRBM for which the marginal distribution over the visible variables is identical to the joint distribution of the initial Markov network.

To give an example of the mapping between Markov networks and sRBMs, consider the Markov network structure in Figure 1(a) and the corresponding sRBM structure in Figure 1(b). For each edge in the Markov network, a hidden variable is added to the sRBM, connecting the visible variables indirectly. Moreover, for some interaction matrix

[TABLE]

where $\theta_{12},\theta_{23}>0$, we have

[TABLE]

Note that the off-diagonal elements in $A$ and $W^{\top}W$ are identical, which is in line with equation (3).

Clearly, the visible variables are mutually conditionally independent given the hidden variables, and vice versa, allowing for parallel sampling of individual variables. Similar to [11], the conditional distributions of the hidden variables given the visible variables are normal:

[TABLE]

Moreover, the conditional distributions of the visible variables given the hidden variables are given by

[TABLE]

where $\sigma()$ is the logistic sigmoid function and $W_{(:,v)}$ is the $v$:th column of $W$. Note that $W_{(vv^{\prime},:)}$ contains only two non-zero elements and that $W_{(:,v)}$ contains as many non-zero elements as there are neighbors of $v$ in $G$.

The mapping between Markov networks and sRBMs provides a straightforward technique for sampling from the joint distribution using a block-Gibbs sampler, which alternates between sampling the visible and hidden variables, and ultimately saving the configurations over the visible variables. Compared to a sequential Gibbs sampler, RBM samplers have been shown to mix slightly slower [11, 12] and they also require an additional step to sample the hidden variables. Still, from a computational perspective, the parallelization-friendly structure of an RBM sampler makes it a very attractive candidate for ultra-high-dimensional Markov network sampling using GPUs or other massively parallel environments.

2.2 Structure learning

Structure learning refers to the process of inferring the graph structure of a model from a set of data generated from that model. The data is assumed to contain $n$ complete i.i.d. joint observations which are represented by an $n$-by-$d$ data matrix $\mathbf{x}=(x^{(1)},\ldots,x^{(n)})$. One of the main challenges of structure learning is to distinguish between direct interactions (or edges) and indirect interactions mediated solely by multi-step paths over other nodes.

Structure learning methods have traditionally been divided into two categories; constraint-based and score-based. Constraint-based methods aim at recovering the model structure using a series of independence tests. Score-based methods formulate the structure learning problem as an optimization problem, which requires a score-function for evaluating how well a graph fits the data and an algorithm for optimizing the score function.

3 Structure learning methods

The methods considered in this work have been developed in various fields spanning statistics, statistical mechanics and computational biology. The methods are from two conceptually very different algorithmic families: pseudo-likelihood methods, which are score-based, and pairwise methods, which are constraint-based using pairwise tests. Common to these methods is that they are applicable to high-dimensional problems and that they provide a way to rank edges by interaction strength. In this section we provide a brief description of the selected methods. For more details, we refer the reader to the original works.

3.1 Pseudo-likelihood methods

The (log-)likelihood function is probably the most important and well-known function in score-based model learning:

[TABLE]

In terms of structure learning, the underlying rationale is that $\theta_{vv^{\prime}}=0$ iff $\{v,v^{\prime}\}\not\in E$. However, maximizing the likelihood of non-chordal Markov networks is computationally intractable, since computing the partition function $Z(\boldsymbol{\theta})$ requires summing over all joint outcomes. Consequently, likelihood-based techniques are limited to small $d$ problems.

As a computationally more convenient alternative, Besag [13] introduced the pseudo-likelihood, which is an approximation of the likelihood consisting of node-wise conditional likelihoods:

[TABLE]

where $-v=V\setminus v$ and $\boldsymbol{\theta}_{v}=(\theta_{v1},\ldots,\theta_{vd})$ is the parameter vector associated with $v$ such that $\theta_{vv}=\theta_{v}$. Initially, optimizing the pseudo-likelihood may seem like a cumbersome task. However, the main computational advantage, compared to the likelihood, is that the intractable partition function $Z(\boldsymbol{\theta})$ cancels out in the conditional distributions. More specifically, the conditional distributions of a node $v$ is given by

[TABLE]

where

[TABLE]

is a local normalizing constant depending on $\boldsymbol{\theta}_{v}$. In addition to being computationally tractable, the pseudo-likelihood function is concave and hence has no local maxima. Also, the maximum pseudo-likelihood estimator is consistent in the sense that the parameter values will approach the true values as $n\to\infty$ [14].

Moreover, examining (4) closer, we see that the conditional distribution takes the form of a logistic regression. From now on, assume that $\theta_{vv^{\prime}}$ and $\theta_{v^{\prime}v}$ represent distinct parameters. Under the equality constraint $\theta_{vv^{\prime}}=\theta_{v^{\prime}v}$, maximizing the pseudo-likelihood can then be thought of as maximizing a sum of coupled logistic regression problems:

[TABLE]

Following [5], we refer to this technique as the symmetric pseudo-likelihood approach. Symmetric pseudo-likelihood methods for structure learning in Markov networks have been developed in [2, 4].

As $d$ grows, optimizing the symmetric pseudo-likelihood eventually becomes too demanding. However, using the above formulation, it is straightforward to apply an even more efficient approach, obtained by simply removing the parameter equality constraint. This results in $d$ decoupled logistic regression problems:

[TABLE]

This approach, which is referred to as asymmetric pseudo-likelihood or regression-based structure learning, is very convenient from a computational perspective, since the subproblems can be solved independently of each other on multiple cores. During optimization, the asymmetric approach considers $\theta_{vv^{\prime}}$ and $\theta_{v^{\prime}v}$ as distinct parameters and will thus (in practice) give out two different estimates for each interaction parameter. Consequently, the solution requires some form of post-processing to be consistent with a Markov network.

In order to avoid overfitting and ill-posedness when maximizing the pseudo-likelihood, it is in general necessary to impose some form of regularization on the model parameters. Next, we will present two asymmetric pseudo-likelihood methods, which are characterized by their type of regularization.

3.1.1 Asymmetric pseudo-likelihood with $L_{1}$-regularization (plmL1)

$L_{1}$-regularization is a popular technique for estimating sparse models from data, since the $L_{1}$-penalty will drive non-informative parameters to zero. Consequently, the $L_{1}$-norm has been considered by several authors as a means to perform structure learning of Markov networks [15, 2, 3]. In particular, [3] introduced the popular asymmetric pseudo-likelihood with $L_{1}$-regularization:

[TABLE]

for which statistical guarantees for consistent graph estimation were proven under certain coherence conditions imposed on the Fisher information matrix. The method is a member of the well-known family of Lasso algorithms, a technique that was pioneered by [16] for covariance selection for Gaussian graphical models. Although the $L_{1}$-norm is not differentiable, the problem is still convex and there exist a wide variety of methods for solving high-dimensional $L_{1}$-regularization problems.

Since $L_{1}$-regularization will actually drive parameters to zero, it is straightforward to infer the graph from the estimated parameters, since a non-zero parameter implies existence of an edge. However, we may obviously come across situations where $|\theta_{vv^{\prime}}|>0$ while $|\theta_{v^{\prime}v}|=0$. The common technique for dealing with the asymmetry issue is to add an edge $\{v,v^{\prime}\}$ if

[TABLE]

or alternatively, if

[TABLE]

The above criteria, which are also known as the $\wedge$ and $\vee$ rules [16], have shown similar accuracy in previous numerical experiments [2, 17].

The strength of the sparsity-inducing regularization term is adjusted by tuning the penalty weight $\lambda$, such that the degree of sparsity increases with a larger penalty weight. By moving across an appropriately chosen range of decreasing $\lambda$-values, one obtains a so-called regularization path, along which the number of non-zero interaction parameters increases from [math] to $d-1$. In other words, the regularization path offers a straightforward way for ranking the edges using the penalty weight.

3.1.2 Asymmetric pseudo-likelihood with $L_{2}$-regularization (plmDCA)

$L_{2}$-regularization is one of the most common techniques for dealing with ill-posed problems in statistics and machine learning, with examples being linear regression [18] and classification with logistic regression or support vector machines [19]. Ekeberg et al. [4] introduced the idea of combining the symmetric pseudo-likelihood for categorical data with $L_{2}$-regularization. The original method was later replaced with an asymmetric version [5]:

[TABLE]

which was shown to give similar results as the original symmetric version. The method was named plmDCA, since it was originally developed for direct coupling analysis (DCA), which is the task of identifying direct couplings, or contacts, between amino-acid positions in protein chains from multiple sequence alignments. PlmDCA is currently one of the leading methods for performing DCA [20]. Being differentiable, the $L_{2}$-norm is a convenient penalty function from the optimization point of view.

In contrast to the $L_{1}$-penalty, $L_{2}$-regularization does not force parameters to zero. Therefore, the regression parameters are turned into an interaction score. The interaction strength of an edge $\{v,v^{\prime}\}$ is calculated as follows. First, the estimated interaction parameters involving $v$ and $v^{\prime}$ are collected into two coupling matrices; one is the result from regressing $v$ on $-v$ and the other from regressing $v^{\prime}$ on $-v^{\prime}$. After adjusting the two coupling matrices by the Ising gauge, they are averaged into a single coupling matrix, of which the Frobenius norm is calculated to give an edge-specific interaction score. Finally, the interaction scores are adjusted by an average-product correction step, which was originally introduced to suppress certain effects related to multiple sequence alignments.

3.2 Pairwise methods

In computational biology, there are a wide range of applications that can be formulated as structure learning problems, one of the most well-known being gene expression networks. In general, the data sets are high-dimensional and of type $n<d$. To tackle the high-dimensional aspect of such data sets, an array of more or less heuristic methods have been developed. Although often lacking any theoretical guarantees that would hold in a general setting, methods of this kind have led to several successful applications.

In this work, we will consider methods that are based on a simple pairwise score, which measures the marginal dependence between two variables. As our pairwise score, we will use mutual information, which is one of the most well-known measures of marginal (or mutual) dependence. The mutual information (MI) between two discrete variables, indexed by $v$ and $v^{\prime}$, is defined by

[TABLE]

In practice, we do not have access to the actual joint distribution $p(X_{v},X_{v^{\prime}})$, which has to be estimated from data. Rather than using the maximum likelihood estimate, we use a smoothing procedure where we add a pseudo-count of $1/4$ to each joint outcome, giving a total pseudo-count of $1$. Thus, our estimates are calculated according to

[TABLE]

where $n(x_{v},x_{v^{\prime}})$ represents the count of the corresponding outcome. From a Bayesian perspective, the estimator can be interpreted as the posterior mean, using a symmetric Dirichlet prior with hyperparameter $1/4$.

The most basic MI-based approach is to the calculate the mutual information for each pair of variables and then rank the strength of the edges according to their MI values (cf. Relevance networks [7]). There is an inherent flaw with this approach, since a pairwise test on its own cannot distinguish between direct and indirect dependencies. Consequently, a variety of post-processing techniques have been designed for dealing with this particular issue.

3.2.1 Context likelihood of relatedness (CLR)

The CLR algorithm was introduced as an approach to remove false positives in gene regulatory networks through a background correction step [9]. After computing the MI between all pairs of nodes, CLR calculates an adjusted interaction score that takes the background distribution of MI values into account. The background distribution of $MI(v,v^{\prime})$ is approximated as a bivariate normal distribution estimated from the MI values involving the nodes $v$ and $v^{\prime}$. More specifically, let $z_{v}(v,v^{\prime})$ and $z_{v^{\prime}}(v,v^{\prime})$ be the z-scores of $MI(v,v^{\prime})$ computed with respect to the MI values involving $v$ and $v^{\prime}$, respectively. The background adjusted interaction score between $v$ and $v^{\prime}$ is then calculated from the z-values according to

[TABLE]

Before calculating the adjusted interaction score, negative z-values are set to zero. As a result, most weight will be given to interactions that stand significantly above the marginal background distributions of the considered variables.

3.2.2 An algorithm for the reconstruction of gene regulatory networks (ARACNE)

ARACNE was proposed as a means for disentangling direct from indirect dependencies when learning gene regulatory networks [8]. The method is based on a notion referred to as the data processing inequality (DPI). The DPI states that if two nodes $v$ and $v^{\prime}$ only interact indirectly through a third node $v^{\prime\prime}$, then

[TABLE]

In other words, information between $v$ and $v^{\prime}$ is in general lost (but never gained) when transmitted through $v^{\prime\prime}$.

ARACNE starts from a graph containing an edge for each non-zero MI value. The algorithm then proceeds by examining each triplet of mutually connected nodes. For each such triplet, the edge associated with the lowest MI value is removed if the difference to the second largest MI value in the triplet is larger than a user-specified tolerance threshold. The triplets are examined independently of each other, meaning that an edge that has been removed with respect to one triplet is still taken into consideration when examining another triplet containing that particular edge, see Figure 2 for an example. The output of the algorithm is therefore not affected by the order in which the triplets are examined.

ARACNE is better suited for certain types of graph topologies. In particular, for a Markov network with a tree-based dependence structure, ARACNE will recover the correct graph for a large enough sample. On the other hand, using a zero tolerance, loops over three nodes cannot be discovered, since the weakest edge in such loops will always be removed (except in the degenerate case). In general, ARACNE is well-suited for networks with locally tree-like structures in which the shortest path dominates the information exchange between two nodes.

3.2.3 Network deconvolution (ND)

The ND algorithm was introduced as a general purpose algorithm for filtering out the indirect effects from an observed correlation matrix containing both direct and indirect effects [6]. By formulating the problem as the inverse of network deconvolution, the algorithm attempts to remove the indirect effects by using eigendecomposition and a result for geometric infinite series of matrices. The fundamental assumption behind the algorithm is that indirect dependencies can be approximated by products of direct dependencies such that the observed marginal dependencies are the sum of direct and indirect dependencies. Although the assumption does not hold in general, the method was shown to perform very well on a collection of real-world problems [6].

Let $G_{\text{dir}}$ and $G_{\text{indir}}$ be symmetric matrices containing the direct and indirect effects of a Markov network. Also, let $G_{\text{obs}}$ be a matrix of observed marginal dependencies. Then ND assumes that

[TABLE]

The goal is to recover the sparse $G_{\text{dir}}$ matrix from the observed $G_{\text{obs}}$. Under the assumption that the maximum absolute eigenvalue of $G_{\text{dir}}$ is strictly smaller than one, the infinite sum can be expressed as $G_{\text{obs}}=G_{\text{dir}}(I+G_{\text{dir}})^{-1}$, which gives the closed form solution

[TABLE]

In practice, the observed dependency matrix $G_{\text{obs}}$ is linearly scaled in order to make the absolute eigenvalues of $G_{\text{dir}}$ smaller than one. The above inverse operation is efficiently computed using the eigendecompositions of $G_{\text{obs}}$ and $G_{\text{dir}}$.

4 Numerical experiments

To empirically compare the accuracy of the methods discussed in the previous section, we performed an extensive numerical simulation study. In addition to the presented methods, we included a plain MI-based method as a baseline. The purpose of the simulation was to investigate how the methods perform in practice in a variety of different scenarios, which have not been tailored for any particular method. In particular, we wanted to examine the relationship between the model size $d$ and the sample size $n$. For more details about the implementation of the methods, see A.

4.1 Experimental setup

To test if the methods are particularly ill- or well-suited for a specific type of network structure, we selected three standard network types with different structural characteristics:

Grid network: A fixed four-nearest neighbour grid graph (also known as a lattice graph).

• 2.

Scale-free network: A random scale-free network generated by the Barabási-Albert model [21]. The algorithm starts from an initial connected network over $m_{0}$ nodes and successively adds new nodes such that each new node is connected to $m$ existing nodes. The algorithm uses a preferential attachment mechanism in which the probability of a node receiving an edge during the generation process is proportional to the number of edges already connected to that node.

• 3.

Small-world network: A random small-world network generated by the Watts-Strogatz model [22]. The algorithm starts by constructing a regular ring lattice in which each node is connected to its $k$ closest neighbors. After the initial step, the algorithm iterates through each existing edge and randomly rewires it with a probability $p$.

For each network type, we generated 50 structures over 40, 200 and 1000 nodes (see A for specific model parameters). Various structural properties of the generated networks are listed in Table 1. All of the generated networks are sparse, with the average number of neighbours varying between 3.35 (grid, $d=40$) and 4.20 (scale-free, $d=1000$). In contrast to the highly regular fixed grid networks, the algorithms behind the randomly generated scale-free and small-world networks are designed to capture certain structural properties observed in many real-world networks. More specifically, the scale-free model generates hub nodes through its preferential attachment mechanism, resulting in large maximum neighbourhoods. In contrast, the small-world model produces high levels of clustering, which is reflected by the average clustering coefficient in Table 1.

For each graph structure, the strength of the associated model parameters were drawn randomly according to

[TABLE]

and the sign of a parameter was determined by a $\text{Bernoulli}(0.5)$ trial. From each model, we generated a data set of size 200, 1000 and 5000 using the sRBM sampler. In total, each method was thus applied on 1350 datasets covering $n/d$-ratios from $1/5$ up to $125$.

4.2 Method evaluation

As a basis for evaluating the performance of the methods, we used precision and recall:

[TABLE]

Precision is the fraction of true edges among the included edges, while recall is the fraction of included true edges among the existing true edges.

Since all methods provide a way to rank the edges according to their estimated interaction strength, we obtain a complete range of precision-recall pairs. Consequently, as our first evaluation measure, we consider the complete precision-recall curve summarized by the Area Under the Curve (AUC). The AUC thus provides a measure of the overall accuracy of the methods across the complete range of recall.

As our second accuracy measure, we use the maximum recall reached by a method under a precision of 0.90 (denoted by $RC_{0.90}$), or equivalently, under a false discovery rate of 0.10. In many applications, it is critical to maintain a high precision for the included edges, since each potential edge needs to be verified manually, or in some cases even experimentally. A high value on $RC_{0.90}$ implies that a method is able to recover a high fraction of the true edges while keeping down the number false positives, in this case, only 1 out of 10 among the discovered edges is a false positive.

4.3 Results

The results of the experiments for the grid, scale-free and small-world model are summarized in Figures 3, 4 and 5, respectively. The upper panels in the figures contain the AUC results and the lower panels contain the $RC_{0.90}$ results. For an overview of the average runtime of the methods, see Table 2. We begin by discussing some more general observations before focusing in on the specific methods.

When comparing the overall accuracy on the different network types, the balanced grid network is the easiest to infer, closely followed by the small-world network, while the scale-free network is the most challenging. Other than that, the results are overall quite consistent across the different networks, however, there is one clear exception. ARACNE shows a drop in accuracy on the small-world network (Figure 5) in comparison to the other methods and networks (Figures 3 and 4). The reason for this is that the small-world model is designed to generate clusters containing a lot of loops of length three (as confirmed by the clustering coefficient in Table 1). As explained in Section 3.2.2, ARACNE is not able to recover loops of length three under a zero tolerance threshold making it particularly ill-suited for the small-world network. As suggested in [8], one could use a non-zero tolerance threshold allowing for some three-loops, however, it is not clear how to choose an appropriate non-zero threshold in practice.

As expected, when comparing results across the rows and columns in the boxplot panels in Figures 3–5, we see that the positive effect on accuracy given by a fivefold increase of the sample size $n$ is larger than the corresponding negative effect of a fivefold increase of the model size $d$. For example, the results found on the diagonal in the boxplot panels correspond to different settings for a fixed $n/d$-ratio of 5. When moving from the upper left corner to the lower right corner, the accuracy is increased.

When comparing the pairwise methods for the smallest sample size $n=200$, they all reached an accuracy comparable to plain MI. For larger samples, however, the post-processing step of the pairwise methods significantly improved the accuracy compared to plain MI (except for ARACNE on the small-world network). The need for filtering out indirect dependencies is particularly clear when examining $RC_{0.90}$. Except for the small-world network, we did not observe any major differences in accuracy between the pairwise methods. In terms of AUC, CLR was marginally more accurate than ARACNE and ND. In terms of $RC_{0.90}$, ARACNE (grid and scale-free network) and ND (all networks) were slightly more accurate than CLR.

When comparing the pseudo-likelihood methods, the $L_{2}$-based plmDCA clearly outperformed the $L_{1}$-based plmL1 throughout this experiment in terms of both evaluation criteria. Overall, plmL1 had a surprisingly low accuracy in this experiment and was even outperformed by plain MI for all but the largest samples. Consequently, we focus on plmDCA in the remaining summary of the simulation results.

When comparing plmDCA and the pairwise methods, the results show that plmDCA requires a relatively large sample to gain an advantage over the pairwise methods through its global approach. Overall, plmDCA and the pairwise methods achieved a comparable accuracy in settings where the $n/d$-ratio was 5 (diagonal in boxplot panels). The pairwise methods achieved a higher accuracy in the small-sample settings where the $n/d$-ratio was 1/5 and 1 (below diagonal in boxplot panels), while plmDCA reached a higher accuracy in the large-sample settings where the $n/d$-ratio was 25 and 125 (above diagonal in boxplot panels). Consequently, in terms of accuracy, simple pairwise methods should be preferred for small samples, while plmDCA should be preferred for large samples. Due to the multifactorial nature of the structure learning problem, it is difficult to give a rule of thumb for what constitutes a large sample size in a more general sense. Still, purely based on these experiments, $n$ should be close to five times larger than the number of variables when $d$ is in the range of 40–1000. When scaling up even further, we would expect the required $n/d$-ratio to come down, however, we cannot verify this in practice simply due to computational reasons.

5 Discussion

We have empirically studied various methods for high-dimensional structure learning of binary pairwise Markov networks. We considered a collection of methods from the families of pseudo-likelihood methods, developed in the fields of statistics and statistical mechanics, and pairwise methods, developed in the field of computational biology. Since the methods have been developed in parallel in different fields, they have not (as far as we are aware) been compared against each other in a context similar to the one considered in this work. Conceptually, pseudo-likelihood methods and pairwise methods are very different. Pseudo-likelihood methods take a global approach, considering all variables simultaneously. In contrast, pairwise methods are based on local pairwise tests, which are post-processed in an attempt to adjust for the global aspect of the problem. Whereas the pseudo-likelihood methods enjoy nice theoretical properties such as statistical consistency, the main motivation behind pairwise methods is typically considered to be computational efficiency.

Although consistency is a desirable property in a structure learning method, it does not guarantee superior performance for finite-sized samples. The numerical results in this work showed that various pairwise methods can compete with, and even outperform, the more elaborate pseudo-likelihood methods, not only in terms of speed, but also in terms of accuracy on large synthetic networks when $n\leq d$, which is a common setting in many real-world structure learning problems. As score-based methods are generally thought to be more accurate than constraint-based methods, this finding is not only surprising, but also a valuable insight from a more practical perspective. In particular, as the size of the considered problems keeps increasing, it becomes ever more important to keep down the computational load of the learning algorithms. A recent example of an ultra-high-dimensional structure learning application (with $n<<d$) is genome-wide epistatic interaction discovery [23, 24], where the networks contain tens of thousands or even up to a hundred thousand single nucleotide polymorphisms.

As part of this work, we also presented an sRBM-based Gibbs sampler, similar to the one proposed in [11]. To further scale up future studies to an ultra-high-dimensional setting, efficient sampling from a ground-truth model is a critical step. Although the sRBM sampler requires an additional step to sample the hidden variables, the main advantage of the sampler, in comparison to the standard sequential sampler, comes from the possibility of sampling the individual nodes in parallel. By exploiting massively parallel platforms, this could enable efficient sampling of ultra-high-dimensional sparse Markov networks over hundreds of thousands of variables.

In this work, we sidestepped the problem of determining a significance threshold for the interaction strength. In practical applications, one is typically only interested in edges that are deemed significant by some statistical criterion. Moreover, the computational efficiency of methods such as ARACNE depends on the possibility to remove edges with non-significant MI values before applying the DPI step. This is the reason behind the seemingly high runtimes of ARACNE compared to the other pairwise methods (see Table 2). The perhaps most common general-purpose method for determining a significance threshold is to use permutation tests in which the values obtained for the original data are compared to values obtained for a corresponding permuted data set, where each column has been randomly rearranged in order to make the variables mutually independent. For example, permutation-based procedures have been used to obtain a significance threshold for mutual information in [7, 8]. Alternatives approaches for determining significance in the model-based pseudo-likelihood methods are based on cross-validation, the extended Bayesian information criterion (eBIC) [17] and inverse finite-size scaling [25].

Acknowledgements

This research was supported by the COIN Centre of Excellence, Academy of Finland grant no. 251170 (JP, YX, SP, MP, JC), KAKENHI no. 17H00764 (YK), and ERC grant no. 742158 (JC). The authors wish to acknowledge CSC – IT Center for Science, Finland, for computational resources.