A global approach for learning sparse Ising models

TL;DR

This paper introduces a global method using $l_1$-regularized logistic regression to simultaneously learn the structure and parameters of sparse Ising models, improving over node-wise approaches by leveraging reciprocal node information.

Contribution

The paper presents a novel global estimation approach for sparse Ising models that estimates all edges and parameters simultaneously, unlike existing node-wise methods.

Findings

The proposed method effectively learns sparse Ising models.

Numerical experiments demonstrate the advantage of the global approach.

The method exploits reciprocal node information during estimation.

Abstract

We consider the problem of learning the link parameters as well as the structure of a binary-valued pairwise Markov model. Under sparsity assumption, we propose a method based on $l_1$- regularized logistic regression, which estimate globally the whole set of edges and link parameters. Unlike the more recent methods discussed in literature that learn the edges and the corresponding link parameters one node at a time, in this work we propose a method that learns all the edges and corresponding link parameters simultaneously for all nodes. The idea behind this proposal is to exploit the reciprocal information of the nodes between each other during the estimation process. Numerical experiments highlight the advantage of this technique and confirm the intuition behind it.

Tables

Table 1. Table 1: Example of a sample of size n 𝑛 n from a generic distribution P ​ ( x 1 , x 2 , … , x p ) 𝑃 subscript 𝑥 1 subscript 𝑥 2 … subscript 𝑥 𝑝 P(x_{1},x_{2},\ldots,x_{p})

	$X_1$	$X_2$	…	$X_p$
$x^{\left(1\right)}$	$-1$	$+1$	…	$+1$
$x^{\left(2\right)}$	$+1$	$-1$	…	$-1$
$x^{\left(3\right)}$	$-1$	$-1$	…	$+1$
$x^{\left(4\right)}$	$-1$	$-1$	…	$+1$
$x^{\left(5\right)}$	$-1$	$-1$	…	$-1$
	⋮	⋮	⋮	⋮
$x^{\left(n\right)}$	$-1$	$-1$	…	$+1$

Equations

$$\begin{array}[]{c}\quad f(x_{i},x_{j}|x_{V\backslash\{i,j\}})=f(x_{i}|x_{V\backslash\{i,j\}})f(x_{j}|x_{V\backslash\{i,j\}})\\ \Updownarrow\\ f(x_{i}|x_{V\backslash\{i\}})\mbox{ does not depend on }x_{j}\end{array}$$

$$X_i\perp X_j|X_{V\backslash \{i,j\}}\iff (i,j)\notin E,$$

$$f(x_1,...,x_p)=\prod _{C\in \mathcal{C}}\varphi (x_C),$$

$$P(x_1,\dots ,x_p)=\frac{1}{Z(\theta )}exp\left(\sum _{i=1}^p{\theta }_i{x}_i+\sum _{(i,j)\in E}{\theta }_{ij}{x}_i{x}_j\right)$$

$$P(x_1,\dots ,x_p)=\frac{1}{Z(\theta )}exp\left(\sum _{(i,j)\in E}{\theta }_{ij}{x}_i{x}_j\right).$$

$$P(x)=\frac{1}{Z(\theta )}exp\left(x^t\Theta x/2\right).$$

$$\begin{array}[]{ccl}P_{\theta_{\cdot r}}(x_{r}|x_{V\setminus\{r\}})&=&\frac{P(x_{V})}{P(x_{V\setminus\{r\}})}=\frac{P(x_{1},\ldots,x_{p})}{P(x_{1},\ldots,x_{r-1},1,x_{r+1}\ldots,x_{p})+P(x_{1},\ldots,x_{r-1},-1,x_{r+1}\ldots,x_{p})}\\ \\ &=&\frac{exp\left(\sum_{i\in N(r)}~{}\theta_{ir}x_{i}x_{r}+\sum_{(i,j)\in E_{-r}}~{}\theta_{ij}x_{i}x_{j}\right)}{exp\left(\sum_{i\in N(r)}~{}\theta_{ir}x_{i}+\sum_{(i,j)\in E_{-r}}~{}\theta_{ij}x_{i}x_{j}\right)+exp\left(-\sum_{i\in N(r)}~{}\theta_{ir}x_{i}+\sum_{(i,j)\in E_{-r}}~{}\theta_{ij}x_{i}x_{j}\right)}\\ \\ &=&\frac{exp\left(2x_{r}\sum_{i\in N(r)}~{}\theta_{ir}x_{i}\right)}{exp\left(2x_{r}\sum_{i\in N(r)}~{}\theta_{ir}x_{i}\right)+1}.\end{array}$$

$$P(X_r=1|x_{V\setminus \{r\}})=\frac{exp\left({\sum }_{i\in N(r)}2{\theta }_{ir}{x}_i\right)}{exp\left({\sum }_{i\in N(r)}2{\theta }_{ir}{x}_i\right)+1}=logistic\left(\sum _{i\in N(r)}2{\theta }_{ir}{x}_i\right),$$

$$\ell ({\theta }_{\cdot r}):=\sum _{i=1}^{n}log{P}_{\theta }(x_r^{(i)}|x_{V\setminus \{r\}}^{(i)})$$

$$P_{\theta }(x_r|x_{V\setminus \{r\}})=\frac{exp\left(x_r{\sum }_{l\ne r}2{\theta }_{lr}{x}_l\right)}{exp\left(x_r{\sum }_{l\ne r}2{\theta }_{lr}{x}_l\right)+1}.$$

$${\hat{\theta }}_{\cdot r}=\arg \underset{\theta \in {\mathcal{R}}^{p-1}}{\min }-\frac{1}{n}\sum _{i=1}^{n}log{P}_{\theta }(x_r^{(i)}|x_{V\setminus \{r\}}^{(i)})+\lambda \Vert \theta {\Vert }_1;$$

$$Q_r^{\star }:={\mathbb{E}}_{{\theta }_{\cdot r}^{\star }}\left\{{\nabla }^{2}log(P_{\theta }(X_r|X_{V\setminus \{r\}})\right\},$$

$${\Lambda }_{min}(Q_1^{\star })\ge C_{min},$$

$$\Vert Q_2^{\star }(Q_1^{\star }{)}^{-1}{\Vert }_{\infty }\le 1-\alpha .$$

$${\theta }_{min}^{\star }:=mi{n}_{i\ne r}|{\theta }_{ir}^{\star }|\ge \frac{10}{C_{min}}\sqrt{d_r}\lambda$$

$${\hat{N}}_{\pm }(r)=N_{\pm }(r),$$

$${\hat{N}}_{\pm }(r):=\{sign({\hat{\theta }}_{ir})i:i\in V\setminus \{r\},{\hat{\theta }}_{ir}\ne 0\}$$

$$N_{\pm }(r):=\{sign({\theta }_{ir}^{\star })i:i\in V\setminus \{r\},{\theta }_{ir}^{\star }\ne 0\}.$$

$$\hat{\theta}^{N-L-m}_{ij}=\hat{\theta}^{N-L-m}_{ji}=\left\{\begin{array}[]{cl}\hat{\theta}_{ij}&if|\hat{\theta}_{ij}|>|\hat{\theta}_{ij}|\\ \hat{\theta}_{ji}&if|\hat{\theta}_{ij}|\leq|\hat{\theta}_{ij}|\end{array}\right.,$$

$$\hat{\theta}^{N-L-M}_{ij}=\hat{\theta}^{N-L-M}_{ji}=\left\{\begin{array}[]{cl}\hat{\theta}_{ij}&if|\hat{\theta}_{ij}|<|\hat{\theta}_{ij}|\\ \hat{\theta}_{ji}&if|\hat{\theta}_{ij}|\geq|\hat{\theta}_{ij}|\end{array}\right..$$

$$\left\{\begin{array}[]{c}\arg\!\min-\frac{1}{n}\sum_{i=1}^{n}logP_{\theta_{\cdot 1}}(x^{(i)}_{1}|x^{(i)}_{V\setminus\{1\}})+\lambda_{1}\sum_{i\neq 1}|\theta_{i1}|\\ \\ \arg\!\min-\frac{1}{n}\sum_{i=1}^{n}logP_{\theta_{\cdot 2}}(x^{(i)}_{2}|x^{(i)}_{V\setminus\{2\}})+\lambda_{2}\sum_{i\neq 2}|\theta_{i2}|\\ \vdots\\ \arg\!\min-\frac{1}{n}\sum_{i=1}^{n}logP_{\theta_{\cdot p}}(x^{(i)}_{p}|x^{(i)}_{V\setminus\{p\}})+\lambda_{p}\sum_{i\neq p}|\theta_{ip}|\\ \end{array}\right.$$

$$\hat{\theta }=\arg \min -\frac{1}{np}\sum _{r=1}^p\sum _{i=1}^{n}log{P}_{{\theta }_{\cdot r}}(x_r^{(i)}|x_{V\setminus \{r\}}^{(i)})+\lambda \sum _{i<j}|{\theta }_{ij}|$$

$$Q^{\star }:={\mathbb{E}}_{{\theta }^{\star }}\left\{{\nabla }^{2}log(P_{\theta }(X)\right\}.$$

$${\Lambda }_{min}(Q_1^{\star })\ge C_{min}$$

$$\Vert Q_2^{\star }(Q_1^{\star }{)}^{-1}{\Vert }_{\infty }\le 1-\alpha .$$

$${\theta }_{min}^{\star }:=mi{n}_{i<j}|{\theta }_{ij}^{\star }|\ge \frac{10}{C_{min}}\sqrt{d}\lambda .$$

$$\hat{S}:=\{ij:{\hat{\theta }}_{ij}\ne 0\}.$$

$$\begin{array}[]{c}P(X_{1}=1|x_{2},x_{3})={exp(\theta_{12}x_{2}+\theta_{13}x_{3})}/{(exp(\theta_{12}x_{2}+\theta_{13}x_{3})+1)}\\ P(X_{2}=1|x_{2},x_{3})={exp(\theta_{12}x_{1}+\theta_{23}x_{3})}/{(exp(\theta_{12}x_{1}+\theta_{23}x_{3})+1)}\\ P(X_{3}=1|x_{1},x_{2})={exp(\theta_{13}x_{1}+\theta_{23}x_{2})}/{(exp(\theta_{13}x_{1}+\theta_{23}x_{2})+1)}.\\ \end{array}$$

$$\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern 0.0pt\kern 2.5pt\kern-5.0pt\left[\kern 0.0pt\kern-2.5pt\kern-5.55557pt\vbox{\kern-0.86108pt\vbox{\vbox{ \halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep& \kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&& \kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr 5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle y$\hfil\kern 5.0pt\crcr\kern 0.86108pt\cr 5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\mbox{\Large$\chi$}_{\cdot 1}$\hfil\kern 5.0pt\cr 5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\mbox{\Large$\chi$}_{\cdot 2}$\hfil\kern 5.0pt\cr 5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\mbox{\Large$\chi$}_{\cdot 3}$\hfil\kern 5.0pt\cr}}}}\right]$}}\quad\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern 0.0pt\kern 2.5pt\kern-5.0pt\left[\kern 0.0pt\kern-2.5pt\kern-5.55557pt\vbox{\kern-0.86108pt\vbox{\vbox{ \halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep& \kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&& \kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr 5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\tilde{X}_{1}$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\tilde{X}_{2}$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\tilde{X}_{3}$\hfil\kern 5.0pt\crcr\kern 0.86108pt\cr 5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\mbox{\Large$\chi$}_{\cdot 2}$\hfil\kern 5.0pt&5.0pt\hfil$\mbox{\Large$\chi$}_{\cdot 3}$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr 5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\mbox{\Large$\chi$}_{\cdot 1}$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$\mbox{\Large$\chi$}_{\cdot 3}$\hfil\kern 5.0pt\cr 5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$\mbox{\Large$\chi$}_{\cdot 1}$\hfil\kern 5.0pt&5.0pt\hfil$\mbox{\Large$\chi$}_{\cdot 2}$\hfil\kern 5.0pt\cr}}}}\right]$}}\quad\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern 0.0pt\kern 2.5pt\kern-5.0pt\left[\kern 0.0pt\kern-2.5pt\kern-5.55557pt\vbox{\kern-0.86108pt\vbox{\vbox{ \halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep& \kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&& \kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr 5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\theta$\hfil\kern 5.0pt\crcr\kern 0.86108pt\cr 5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\theta_{12}$\hfil\kern 5.0pt\cr 5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\theta_{13}$\hfil\kern 5.0pt\cr 5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\theta_{23}$\hfil\kern 5.0pt\cr}}}}\right]$}},$$

$$\begin{array}[]{cl}P(Y=1|\tilde{x_{1}},\tilde{x_{2}},\tilde{x_{3}})&=logistic(\theta_{12}\tilde{x_{1}}+\theta_{13}\tilde{x_{2}}+\theta_{23}\tilde{x_{3}})\\ &=exp(\theta_{12}\tilde{x_{1}}+\theta_{13}\tilde{x_{2}}+\theta_{23}\tilde{x_{3}})/(exp(\theta_{12}\tilde{x_{1}}+\theta_{13}\tilde{x_{2}}+\theta_{23}\tilde{x_{3}})+1)\end{array}.$$

$$Accuracy=(TP+TN)/(TP+TN+FN+FP),$$

Full text

1 Introduction

Ising models are fundamental undirected binary graphical models that capture pairwise dependencies between the input variables. They are well-studied in the literature and have applications in a large number of areas such as physics, computer vision and statistics ([7]; [10]; [14]). One of the core problems in understanding graphical models is structure learning, that is, recovering the structure of the underlying graph given a random sample from the population distribution. The literature is very rich and it is difficult to make a complete picture of this topic. In the seminal paper, [4] proposed a greedy algorithm for undirected graphical models assuming the underlying graph is a tree. Since then, various algorithms have been proposed for structure estimation of graphical models. They can be broadly classified into two categories: combinatorial algorithms and convex relaxation algorithms. The first category of methods are based on an empirical estimate of some statistics by which performing local tests on small groups of data, and then combining them to output a graph structure: among many others [2] and [3] ultimately propose an algorithm based on the empirical frequencies, [1] propose an algorithm based on empirical conditional variation distance, and [9] proposes a sparsitron algorithm based on a weighted voting scheme. All these methods focus on graph estimation but not on parameters estimation. On the other hand, the second category includes methods based on penalized convex optimization. Among many others, under different assumptions on the structure of the underlying graph [12], [15] and [13]. In particular, [13] was a pioneering paper in this direction proposing to learn the graphical structure by a node-wise approach using an $\ell_{1}$-penalized logistic regression. In that paper for the first time the authors were also interested in estimating the link parameters, or more precisely their sign. At the same time, in [8], exploiting the same idea, attention was paid to both aspects, namely the estimation of the graph structure as well as the estimation of the values of the link parameters and not only their sign.

We start from these last two works and propose a modification of the method in [13], that allows us to estimate the structure of the graph and the link parameters (not only their sign) in a global manner and not one node at a time. The proposed global approach takes the form of a parallelization of the node-wise approach, in the sense that it learns all the nodes simultaneously, thus exploiting the mutual information of each node on the others.

The paper is organized as follows. In Section 2 we present the mathematical and statistical formulation of the problem. In Section 3 we present the node-wise approach which represents the starting point of the proposed method, and in the same section we describe our global approach. In the last Section, we show a simulation study comparing the node-wise and the global approach confirming the intuition behind it.

2 Mathematical framework: the Ising model

For a complete and exhaustive treatment of graphs theory we refer to [11]; below we give only definitions and properties necessary for this work. A finite graph $G=(V,E)$ consists of a finite collection of nodes $V=\{1,2...,p\}$ and a collection of edges $E\subseteq V\times V$. For the scope of this work, we consider graphs that are undirected, namely graphs whose edges are not ordered, i.e. there is no distinction between the edges $(i,j)$ and $(j,i)\in E$. Moreover, for any $i\in V$, $N(i):=\{j\in V:(i,j)\in E\}$ is the set of neighbours of node $i$ and $C\subset V$ is a clique if $(i,j)\in E$ for all $i,j\in C$ such that $i\neq j$.

In this paper the notion of a graph is used to keep track of the conditional dependence relationship between random variables of a complex system. By complex system here we mean a jointly distributed random variables $(X_{1},X_{2},...,X_{p})$ that interact with each other. In the following, a formal definition of conditional independence relationship is given:

Definition: Two random variables ($X_{i}$, $X_{j}$) of a random vector $(X_{1},\ldots,X_{p})$ are conditionally independent, $X_{i}\perp X_{j}|X_{V\backslash\{i,j\}}$, if

[TABLE]

where $f(\cdot)$ stands for density distribution or probability mass function and $x_{S}:=(x_{s},s\in S)$, with $S$ any subspace of $V$.

Associated with an undirected graph $G=(V,E)$ and a system of random variables $X_{V}$ indexed in the vertexes set $V$ there is a range of different Markov properties which establish how much the graph is explanatory of the conditional independence property of the random variables, see [11]. Specifically, in this work we deal with system of random variables which are pairwise Markov with respect to an undirected graph $G=(V,E)$, so it holds that

[TABLE]

which establishes conditional independence among two variables $X_{i}$ and $X_{j}$ if and only if their corresponding nodes in the graph $G$ are not connected. Another way an undirected graphical model can encode the conditional dependency relationships between the system variables is through the factorization property. Let $\mathcal{C}$ be the set of all possible cliques in a graph $G$, then the distribution factorizes as

[TABLE]

where $\phi(X_{C})$ is a positive function (potential function) which depends only on the variables corresponding to the nodes in $C$, $\forall C\in\mathcal{C}$. Under the hypothesis that the joint distribution $f(x_{1},...,x_{p})$ is positive, the pairwise Markov property and the factorization property are equivalent as claimed by the Hammersley-Clifford theorem proved in ([11] cfr. Theorem 3.9).

In this paper, we work under such hypothesis, in particular with a system of binary random variables with values in $\{-1,1\}$ and for which the multivariate distribution (probability mass function in such a case) factorizes in the following way:

[TABLE]

with $\theta_{i}$ and $\theta_{ij}\in\mathbb{R}$ some parameters and $Z(\theta)$ a constant making the probabilities sum to one (usually called partition function), see [11]. The set $E$ represents the edges’ set of the undirected graph which is pairwise Markov with respect to the distribution and it represents the conditional dependency relationships among the system variables. This model is known as Ising model and it is used in many application of spatial statistics such as modelling the behaviour of ferromagnets, since in such case the discrete variables represent magnetic dipole moments of atomic spins arranged in a graph that can be in one of two state $\{+1,-1\}$, see [13]. In particular, in this paper we consider model with no first order terms, which corresponds to have a zero external field interacting with the ferromagnets’ system:

[TABLE]

From a probabilistic point of view, having no first order terms makes the model symmetric to switching the value of the variable in all the graph nodes, i.e. $P(x_{1},\ldots,x_{p})=P(-x_{1},\ldots,-x_{p})$. This is the model considered in [13].

Note that, since in graph $G=(V,E)$ there is no distinction between edge $(i,j)$ and $(j,i)$, in Equation (4) there is no distinction between parameter $\theta_{ij}$ and $\theta_{ji}$, from a physical point of view this mean that the interaction strength between two variables is a reciprocal/symmetric relationship. For convenience of exposition, we define the vector $x=(x_{1},\ldots,x_{p})$ and the symmetric zero diagonal matrix $\Theta$ by which it is possible to express the joint probability in (4) by the following formula

[TABLE]

Matrix $\Theta$ is our unknown parameter; its support encoding the edges’ set $E$ of the graph, its $p(p-1)/2$ upper off-diagonal elements representing the link parameters.

Our perspective is inferential. Given a sample of size $n$ extracted from an unknown distribution $P(x)$ of form (4), as the one represented in Table 1, we want to learn both the structure of the undirected graph $G$ which is Markov with respect to the distribution as well as the link parameters itself. In other words, we want to learn the dependence/independence conditional relationships between the system variables (i.e. the pairwise Markov graph structure) as well as the strength of these relationships (i.e. the numerical value of link parameters $\theta_{ij}$). Both of these objectives are achieved if matrix $\Theta$ is correctly estimated, i.e. if the whole vector $(\theta_{ij})_{i<j}$ of its upper off-diagonal elements is correctly estimated.

3 Estimation methods

In the milestone paper [13] the authors propose one of the most efficient methods for recovering the graph $G$ associated to an Ising model as well as the sign of its link parameters, $sign(\theta_{ij})$. Their method does not make hypotheses on the graph’s structure apart from the sparsity. From that work many others have emerged in the literature, lastly [16]. In the following we revise it to make our proposal clearer.

Let us consider the conditional distribution of one system variable with respect to the others. Let us fix a variable $X_{r}$, with $r\in\{1,2,..p\}$ and consider the set of the remaining variables $X_{V\setminus\{r\}}$; define $E_{-r}=E\setminus\{(i,r):i\in N(r)\}$ as the subset of edges not involving the $r$-th node, hence we have

[TABLE]

Evaluating the previous formula for $x_{r}=1$, it follows that

[TABLE]

where $logistic(*)=e^{*}/(e^{*}+1)$. It is then possible to learn vector $\theta_{\cdot r}$ solving a $l_{1}$-penalized logistic regression problem, where $X_{r}$ and $X_{V\setminus\{r\}}$ play the role of response variable and covariates respectively.

In particular, let $\{x^{(i)}\in\{+1,-1\}^{p}\}_{i=1,...,n}$ be the sample of size $n$ extracted from the unknown population, then the loglikelihood function for the unknown parameter vector $\theta_{\cdot r}=(\theta_{lr})_{l\neq r}$ is the following function

[TABLE]

with

[TABLE]

Since we expect that the elements of the unknown vector corresponding to the neighbours of node $r$ are different from zero (i.e. $\theta_{lr}\neq 0$ $\leftrightarrow$ $l\in N(r)$), while the elements corresponding to the non-neighbours of node $r$ are zero (i.e. $\theta_{lr}=0$ $\leftrightarrow$ $l\notin N(r)$), the following $l_{1}$-penalized logistic regression problem is solved

[TABLE]

$N(r)$ being the true set of active variables into problem in Eq. (6) with cardinality $d_{r}=card(N(r))$ equals the degree of node $r$. Let us revise the theoretical property of this estimator, as exposed in the original paper [13]. Denote the true model parameter $\theta^{\star}$. For any fixed node $r\in V$, define the Hessian matrix ($p-1\times p-1$) of the likelihood function $\mathbb{E}(log(P_{\theta}(X_{r}|X_{V\setminus\{r\}}))$ as evaluated at the true model parameter $\theta^{\star}_{\cdot r}$

[TABLE]

and denote $Q^{\star}_{1}$ and $Q^{\star}_{2}$ the restriction of $Q^{\star}_{r}$ to the index sets $N(r)\times N(r)$ and $V\backslash\{N(r)\cup\{r\}\}\times N(r)$ respectively, i.e. $Q^{\star}_{1}$ and $Q^{\star}_{2}$ have dimension $d_{r}\times d_{r}$ and $p-d_{r}-1\times d_{r}$, respectively. We can now state the following three assumptions:

1. A1)

Dependency condition: there exists a constant $C_{min}$ such that

[TABLE]

with $\Lambda_{min}$ the minimum eigenvalue. 2. A2)

Incoherence condition: there exist $\alpha\in(0,1]$ such that

[TABLE] 3. A3)

Minimum strength signal:

[TABLE]

Under assumptions A1), A2) and A3), Theorem 1 in [13] proves that, when $\lambda=\mathcal{O}(\sqrt{logp/n})$ and $n$ is big enough, with high probability the solution $\hat{\theta}_{\cdot r}$ of problem (6) consistently estimates the signed set of node $r$ neighborhood, i.e.

[TABLE]

with

[TABLE]

The reason why the authors do not consider $\hat{\theta}_{ir}$ as estimate of the link parameter lies in the fact that they solve problem (6) for each node $r\in V$ independently of the other nodes so that $\hat{\theta}_{ir}\neq\hat{\theta}_{ri}$. This means that, if for each variable $X_{r}$ we place the estimate of its link parameter vector $\hat{\theta}_{\cdot r}$ into the extra diagonal elements of the corresponding column of matrix $\hat{\Theta}$, then this leads to a non-symmetric estimated matrix. For this reason in [8] two procedures for symmetrizing this method are proposed. The first procedure, named N-L-m (Neighborhood-based Logistic minimum), works in the following way:

[TABLE]

similarly, the second procedure, named N-L-M (Neighborhood-based Logistic Maximum), works in the following way:

[TABLE]

We stress that, both procedures consist of two steps: in a first step for each node $r\in V$ vector $\hat{\theta}_{\cdot r}$ is estimated by solving (6), in a second step procedure (8) and (9) is applied respectively.

3.1 Global Logistic method

In this paper we propose a new method for learning both the graph $G=(V,E)$ as well as the link parameters $(\theta_{ij})_{i<j}$. The idea behind our proposal is to learn matrix $\Theta$ globally instead of column-wise as in the previous two methods. Our idea can be seen as a kind of parallelization of the node-wise approach, in the sense that we learn all the nodes simultaneously, thus using the reciprocal information of each node on the others. In fact, when solving problem (6) for node $r$, parameter $\theta_{ir}$ also come into play when solving problem (6) for node $i$. This fact is more clear if we write the $p$ independent problems in one system of problems:

[TABLE]

We stress that the solution of system (10) are $p$ vectors each of length $p-1$, hence a total of $p(p-1)$ parameters are obtained instead of $p(p-1)/2$ of real interest. Our proposal is to solve the following problem:

[TABLE]

where $P_{\theta_{\cdot r}}(x^{(i)}_{r}|x^{(i)}_{V\setminus\{r\}})$ is given in eq. (5) and the number of unknown parameters is exactly $p(p-1)/2$. The same idea of jointly infer the neighbourhood of all the graph nodes simultaneously has been proposed in [6] for the problem of learning Gaussian Graphical Models. To highlight the fact that this method learns the graph in a global fashion and not node-wise, we call this procedure Global Logistic (G-L). We can now state a theoretical property of this estimator, by exploring results for a single node logistic regression in [13]. Define the Hessian matrix ($p(p-1)/2\times p(p-1)/2$) of the likelihood function $\mathbb{E}(log(P_{\theta}(X)))=\mathbb{E}(\sum_{r=1}^{p}log(P_{\theta_{\cdot r}}(X_{r}|X_{V\setminus\{r\}}))$ as evaluated at the true model parameter $\theta^{\star}:=(\theta^{\star}_{ij})_{i<j}$

[TABLE]

Define the row/column index set $S:=\{ij:(i,j)\in E\}$, note that it equals the set of active variables for problem in Eq. (11), i.e. $S=\{ij:\theta^{\star}_{ij}\neq 0\}$ and, moreover, define $Q^{\star}_{1}$ and $Q^{\star}_{2}$ the restriction of $Q^{\star}$ to $S\times S$ and $S^{c}\times S$ respectively, i.e. $Q^{\star}_{1}$ and $Q^{\star}_{2}$ have dimension $d\times d$ and $p(p-1)/2-d\times d$, respectively with $d=card(S)$. We can now state the following three assumptions, similar to those in the previous section :

1. B1)

Dependency condition: there exists a constants $C_{min}$ such that

[TABLE] 2. B2)

Incoherence condition: there exist $\alpha\in(0,1]$ such that

[TABLE] 3. B3)

Minimum strength signal:

[TABLE]

Under assumptions B1), B2) and B3), by using Theorem 1 in [13], we can state that when $\lambda=\mathcal{O}(\sqrt{log(p(p-1)/2)/np})$ and $n$ is big enough, with high probability the solution $\hat{\theta}$ of problem (11) consistently estimates the set of active variables, i.e. $\hat{S}=S$ with

[TABLE]

This result ensures the consistency of the method at least in terms of support of the matrix $\Theta$, but its assumptions are neither verifiable nor usable from a practical point of view. It remains however the possibility to use the choice of $\lambda$ suggested by this result, even if not optimal, to make a first estimate of the matrix $\Theta$. This is, in fact, the choice of $\lambda$ made by the same authors in [13] to perform their numerical tests, and therefore, this is the choice that we make to perform our experiments in the next Section. In the following Section, in fact, we show some numerical tests to understand how the proposed procedure improves with respect to the N-L-m and N-L-M presented above. Finally, it is worthwhile to observe that, in principle from computational point of view, there is almost no difference in solving $p$ independent $l_{1}$-penalized logistic regression problems of size $n$ and dimension $p-1$ each, and solving one single $l_{1}$-penalized logistic regression problem with size $np$ and dimension $p(p-1)/2$; however this can be not properly true if a standard software package is used to solve the optimization problem, without exploiting the particular structure of the likelihood function. On the other hand, there is difference in selecting $p$ different regularization parameters ($\lambda_{1},\ldots,\lambda_{p}$) when applying N-L-m and N-L-M procedures, and selecting one single parameter $\lambda$ when applying G-L procedure. Hence, in principle, this is an important advantage of the global approach with respect to the node-wise approach.

4 Numerical experiments

In this Section, we show some numerical experiments to study the performance of the proposed method. Before presenting results it is necessary to specify how we generated data, what algorithm we used, how we chose the regularization parameter and what kind of indexes we used to measure performance.

Data generation: let us denote $\chi$ as data matrix of dimension $n\times p$ which represents $n$ independent realizations of an Ising model of type (4). This means that we fix a matrix $\Theta$ (to be specified later) and we generate samples from distribution given in (4). While in the case $p\leq 10$ this is possible for any computer resource, in the case $p>10$ the set of possible realizations of system variables, $2^{p}$, may became prohibitive. For this reason, we describe the Gibbs sampling scheme we have adopted to generate samples from distribution (4) when its dimension $p$ is high.

randomly choose an initial state:

• $x^{(0)}\in\{+1,1\}^{p}$

for k=1,… generate a new sample $x^{(k)}$ by the following procedure:

• for $r=1:p$
• sample $x_{r}^{(k)}$ from $P(x^{(k)}_{r}|x^{(k)}_{1},...,x^{(k)}_{r-1},x^{(k-1)}_{r+1},...,x^{(k-1)}_{p})$ given in (5)
• end

of course the $n$ samples are taken after an appropriate burn-in time.

Algorithm: we use the Matlab function lassoglm to solve both problem (10) and problem (11), see [5] for details. While for solving problem (6) we use the $r$-th column of data matrix \mbox{\Large\chi}_{\cdot r} as regressor and the remaining columns as covariates, for solving problem (11) the construction of regressor and covariates is more complicated. For completeness we describe it in the simplest case $p=3$ for which we have:

[TABLE]

In this case the regressor vector, the covariates matrix and the unknown parameters vector are the following:

[TABLE]

being $y$ a sample of size $3n$ of a binary random variable $Y$ and $p(p-1)/2=3$ the effective number of unknown link parameters. With this choice indeed we have that

[TABLE]

Choice of $\lambda$: in order to be fair in a comparative study it is legitimate to fix regularization parameter $\lambda=\sqrt{log(p)/n}$, with $p$ number of unknown parameters and $n$ sample size, which is known from the results revised in Sections 3 to be order of the optimal regularization parameter. The same choice is indeed done in the numerical experiments of both [13] and [8] respectively. Hence, in our experiments in each logistic regression of (10) we fix $\lambda=\sqrt{log(p-1)/n}$, while in logistic regression (11) we fix $\lambda=\sqrt{log(p(p-1)/2)/pn}$. Of course, this is not the optimal choice of $\lambda$, which should be performed by a model selection procedure like, CV (Cross Validation), BIC (Bayesian Information Criterion) or MDL (Minimum Description Lenght); however, such a problem is not pursued here due its computational demanding effort.

Performance indexes: since we are interested both in reconstructing the structure of the graph and in estimating the link parameters, we calculate two different indexes of performance. The first index measures how the method correctly estimates the structure of the graph and it is defined as:

[TABLE]

where

• $TP$ is the number of edges present in the graph and correctly identified (i.e. $\theta_{ij}\neq 0\wedge\hat{\theta}_{ij}\neq 0$),
• $TN$ is the number of edges not present in the graph and correctly identified (i.e. $\theta_{ij}=0\wedge\hat{\theta}_{ij}=0$),
• $FN$ is the number of edges present in the graph and not correctly identified (i.e. $\theta_{ij}\neq 0\wedge\hat{\theta}_{ij}=0$) and
• $FP$ is the number of edges not present in the graph and not correctly identified (i.e. $\theta_{ij}=0\wedge\hat{\theta}_{ij}\neq 0$).

In all previous definitions $\hat{\theta}$ is the estimator obtained in Eqs (8), (9) and (11) respectively. Note that the measure in (13) is a scaled measure inherited from the binary classification literature, $0\leq accuracy\leq 1$.

The second index measures how the method correctly estimates the link parameters and it is defined as the $\ell_{2}$-norm of the difference between the true and estimated parameters vector:

[TABLE]

We can now describe the specific settings we chose for numerical experiments. We generated matrix $\Theta$ according to the mixed coupling model, meaning that $\theta^{\star}_{ij}=\pm 0.5$ with equal probability as described in Section 6 of [13]. We propose three examples, meaning that we generated three matrices $\Theta$, of different sizes ($p=25,50,100$) and different sparsity regimes, specifically 12%, 6% and 3% are the percentages of non-zero elements respectively. For each of the three examples we consider four different sample sizes mimicking different dimensional regimes. In Figures 1, 2 and 3 we plot the $accuracy$ obtained on 20 independent runs by all the three methods N-L-m, N-L-M ,G-L for the three examples respectively. We can observe that, as expected, $accuracy$ improves when the sample size increases for all the three methods, however the mutual proportions between the $accuracy$ of the different methods preserve in almost all the cases. It is important to stress that, accuracy is not optimal in such numerical experiments since the choice of $\lambda$ is fixed a priori for all the methods; to improve $accuracy$ a proper model selection technique should be implemented, but for comparison purpose this is not necessary. The same choice of $\lambda$ was made in [13]. We can observe that in terms of accuracy the G-L approach slightly improves with respect to N-L-m and N-L-M, but this is not surprising since the methods work under the same hypothesis and almost with the same estimator. Similar observations can be drown, with respect to the estimation error, i.e. the ability of estimator in reconstructing the value of the link parameters. In Figures 4, 5 and 6 we plot the corresponding $Err$ index obtained on the same independent runs and settings. The advantage of G-L with respect to N-L-m and N-L-M is now more evident across all the examples and across all the sample sizes. The ability of the G-L estimator to better estimate the link parameters is due to the better use of the information, because it focuses on the exact number of unknown parameters. This also explains why, increasing the sample size, this difference decreases: as the number of data increases, information increases so that the advantage of G-L with respect to N-L-m and N-L-M is less evident.

5 Conclusions

We presented a modification of the node-wise logistic regression estimator, which learns the graph jointly and not one node at a time, exploiting the reciprocity of the relationship between nodes, that is the symmetry of the unknown matrix $\Theta$. Within the limits of the numerical experimentations carried out in this work, we can say that this modification slightly improves the ability of the method to reconstruct the graph but it imporves quite significantly the estimation of the link parameters. Both of these improvements are reduced as the sample size increases, decreasing the difficulty of the problem. The proposed technique is based on the hypothesis of sparsity, which cannot be ignored. In the case of no sparsity hypothesis, both node-wise and global approaches should be applied without the Lasso penalty, but necessarily with a number of data $n>p$ that allows the likelihood maximization without constraints.

Possible future developments concern the writing of an ad hoc algorithm for the proposed problem that exploits the particular form of the likelihood to accelerate its convergence. As well as, possible extension to model with interaction terms of order greater than two.

The Matlab codes used to produce results of this paper are available at http://www.iac.cnr.it/ danielad/software.html.

Acknowledgments

This work was supported partially by CNR Italian Flagship project InterOmics and partially by INdAM-GNCS Project 2018.