The problem of graphical selection is to correctly estimate the graph structure of a Markov random field given samples from the underlying distribution. This problem is central to statistics and machine learning, with consequences for a variety of application domains, e.g., image analysis [2], social networks, and computational biology [6].

In general, a graphical model is a structured representation of a joint distribution of potentially dependent random variables. In the models considered in this paper, we use the vertices of a graph G=(V,E) to represent our (binary) random variables. Let |V|=p. Corresponding to the graph, we associate a joint probability distribution on the p random variables with the edges E determining the exact joint distribution. Perhaps a big strength of the graphical model representation is the visualization possible about which models are close to each other. One would expect distributions with similar graphs to be close to each other and indeed, as an auxillary result, we formalize this notion in this paper.

The distributions we consider are Markov random fields with respect to their graphs. Namely, each distribution satisfies a spatial Markov-type independence relationship—any random variable X, conditioned on those adjacent to it in the graph associated with the distribution, is independent of all other variables.

Consider a set G of such models. A sample from a distribution in G is a p-dimensional binary vector, corresponding to the realization of the p random variables. Given several samples from a distribution, formally, the graph selection problem requires us to choose the structure of the model in G that best fits the data.

If the underlying graph is known to be tree-structured, then the problem can be reformulated as a maximum-weight spanning tree problem [5], and solved in polynomial time. On the other hand, for fully general graphs with cycles, the problem is known to be computationally hard [4]. Nonetheless, a variety of practical methods have been proposed, including constraint-based approaches [14], heuristic search, and ℓ1-based relaxations [9], [16], [17]. In particular, it can be proven that some of the ℓ1-based approaches are consistent for model selection under particular scalings of the graph size, degrees, and number of samples [9], [16], and for the case where each degree is bounded by d, the sample complexity has been independently obtained in [3].

Of complementary interest—and the focus of the paper—are the information-theoretic limitations of the graphical selection problems for pairwise binary Markov random fields, also known as Ising models [1]. The Ising model (1) has its origins in statistical physics [1], where it is used to model physical phenomena such as crystal structure and magnetism. Its relative simplicity and the ability to capture practically useful dependencies makes it a candidate for image processing [2], [7], gene network analysis, analysis of social networks, and recently in coding theory starting from [13] (see e.g. [10] for details and a list of references on the topic) and communications, e.g. [8].

At the same time, the dimensionality of the data sets—the number of random variables in question—has increased significantly for many practical problems. For example, gene networks attempt to model interactions between the expression levels of thousands of genes, usually using only tens of samples. The natural question, one that will be addressed here, is whether it is even possible to obtain anything useful from such limited data, and if so, what?

We are given n i.i.d. samples—namely np-dimensional vectors—from a fixed but unknown model in G. As mentioned before, we consider infering E of the graph alone—namely given samples from an unknown model in G, we have to infer the structure (edge set E) of the unknown model. This setting is useful for problems in gene regulatory network inference, for example, where there is a severe shortage of data to do a model selection problem.

How must n scale with parameters corresponding to the graphs of the G so that there exists a method that recovers the correct probability model? Conversely, for what scalings does any method fail to recover the underlying structure correctly? Note that from this perspective, the graphical selection problem is a channel coding problem, in which the messages are the graphs in our family, and each use of the channel corresponds to taking an i.i.d. sample.

We analyze the information-theoretic limitations of this problem under high-dimensional scaling, when the connections of each variable in the underlying MRF is bounded by d. We derive both necessary and sufficient conditions on the scaling of the triplet (n,p,d) for asympotically reliable reocovery of the graph structure.

In Section II, we set up the problem of determining the sample complexity of graphical selection formally. Section III contains the formal statements of our sample complexity results. Section VI contains our results regarding the divergence between Ising models, showing that, roughly speaking, distributions represented by graphs that are similar are indeed closer. Section V outlines the core results that we use for the proofs of results in Section III, and Sections VII and VIII outline the proofs.

We begin with a precise statement of the problem.

Given an undirected graph G=(V,E), let us associate with each vertex i∈V a binary random variable Xi∈{−1,+1}. We then consider a distribution over the random vector X=(X1,…,Xp) with probability mass function

View Source{\BBP}_{\lambda}(x)={1\over Z(\Lambda)}\exp\left\{\sum_{(s, t)\in E}\lambda_{st}x_{s}x_{t}\right\}, \eqno{\hbox{(1)}}

We will refer to the (p2) parameters, λ{st},s,t∈V., as the edge parameters.

At the outset, the sample complexity of inferring an edge set depends in on some form of an upper bound on the edge parameters. As the following example indicates, if the edge parameters are unboundedly large, it may be impossible to distinguish distinct models (potentially with distinct edge sets), even given infinite amounts of data.

To take into account the above observation, our results for a class G depend on B —the smallest real number that satisfies the following for all Λ∈G. If Λ is associated with the graph (V,E), with parameters λst on (s,t)∈E, then for all s∈V,

View Source\sum_{(s, t)\in E}\vert \lambda_{st}\vert \leq B. \eqno{\hbox{(2)}}

No matter what the class G is, selection of a graphical model is not necessarily the same as simply estimating pairs of vertices (s,t) for which cov{Xs,Xt}≠0; in-deed, given the distribution (1), it can be seen that Xs and Xt could be correlated even when there is no direct edge connecting them.

On the other hand, if each edge parameter that is nonzero is finite, given infinite amounts of data, two distinct models Λ(i) and Λ(j) can always be told apart by looking at the vector of all covariances {EXsXt}(s,t)∈V2s≠v.

We observe that a new elementary proof of the above statement can be constructed along the lines of the arguments presented here (see [11] for full proofs). Morever, we note that a more general result on exponential families exists along these lines, see [15].

With finite data covariances can only be approximated. The question therefore is, how fine should these approximations be? This question determines the number of samples that are necessary for the graph selection problem, and a constructive argument determines the number of samples that are sufficient.

In [12], we considered fitting bounded degree models.

Let θ be a (p2) length vector of real numbers. We index the components of the vector by a pair of numbers {s,t},s,t∈V, and write the components as θst. Furthermore, assume that for some real number λ>0,|θ|≥λ1, where the vector inequality is taken to mean a component-wise one. Then the set Gd,p,θ contains all models on vertices V such that

1. each vertex v in V has degree at most d, namely for the edge set E associated with the model, |{e∈E:v∈e}|≤d, for some d≥1

2. if (s,t)∈E, the edge parameter is θst, else 0.

We build on the above case to let the edge parameters be arbitrary, as long as they are bounded. Specifically,

View Source{\cal G}_{d, p}=\cup {\cal G}_{d, p,\theta},

One observation regarding the models described above: in a model belonging to any of the classes above, for all vertices s∈V,EXs=0.

We begin by stating the results obtained for the sample complexity of inferring the structure of graphical models. The following theorem addresses the sufficiency part.

Theorem 1:

Let 0<δ≤1 and let

λ=min(s,t)∈V2λst.

View Source\lambda=\min_{(s, t)\in V^{2}}\lambda_{st}.

If

n>(3e2B+1)2sinh4(λ2)(2logp+log1δ)

View Sourcen>{(3e^{2B}+1)^{2}\over \sinh^{4}({\lambda\over 2})}\left(2\log p+\log{1\over \delta}\right) then ∃ decoder q:Xn→Gd,p such that for all Λ(i),1≤i≤N,

P(q(Xn)≠Λ(i)|Λ(i))<δ.□

View Source{\BBP}(q(X^{n})\neq\Lambda^{(i)}\vert \Lambda^{(i)})<\delta.\qquad\qquad\quad \square

The following theorem addresses the number of samples that are necessary for the estimation of the edge structure.

Theorem 2:

For all decoders q:Xn→Gd,p, if B≥1 and

n≤12eB/2dλ(logpd−1)32(eλ−e−λ)+14dlogpd

View Sourcen\leq{1\over 2}{e^{B/2}d\lambda(\log pd-1)\over 32(e^{\lambda}-e^{-\lambda})}+{1\over 4}d\log{p\over d} there ∃ model Λ(i) such that

P(q(Xn)≠Λ(i)|Λ(i))≥12.

View Source{\BBP}(q(X^{n})\neq\Lambda^{(i)}\vert \Lambda^{(i)})\geq{1\over 2}.

To quantify how far apart and therefore, how easy it is to tell models apart, we define distance measures between models. In Section IV-B, we examine the behavior of the edge inference problem for various values of the parameter B.

A. Properties of the distance measures

We use the following distance between models Λ(i) and Λ(j)

D(Λ(i),Λ(j))=defD(Λ(i)+Λ(j)2∥Λ(i))+D(Λ(i)+Λ(j)2∥Λ(j)),(3)

View Source\eqalignno{{\cal D}(\Lambda^{(i)}, \Lambda^{(j)})&\buildrel{\rm def}\over {=}D\left({\Lambda^{(i)}+\Lambda^{(j)}\over 2}\Vert \Lambda^{(i)}\right)+\cr &\qquad D\left({\Lambda^{(i)}+\Lambda^{(j)}\over 2}\Vert \Lambda^{(j)}\right),&\hbox{(3)} } where the (Λ(1)+Λ(j))/2 denotes the model obtained by averaging the parameters of the two models. This bounds the Chernoff exponent in the hypothesis test between models Λ(i) and Λ(j). In Theorem 3, we bound D(i,j) for any two models.

Another useful measure of distance we consider is the symmetrized KL distance between two models, Λ(i) and Λ(j),

J(Λ(i),Λ(j))=defD(Λ(i)∥Λ(j))+D(Λ(j)∥Λ(i)).(4)

View Source\eqalignno{{\cal J}(\Lambda^{(i)}, \Lambda^{(j)})&\buildrel{\rm def}\over {=}D(\Lambda^{(i)}\Vert \Lambda^{(j)})+\cr &\qquad\quad D(\Lambda^{(j)}\Vert \Lambda^{(i)}). &\hbox{(4)}}

Clearly D(Λ(i),Λ(j))≥0 and J(Λ(i),Λ(j))>0. For the above defined distances between models Λ(i) and Λ(j), we write D(i,j) and J(i,j) where there is no ambiguity. Note that

J(i,j)=J(Λ(i),Λ(j))=∑(s,t)∈Ei∪Ej(λist−λjst)(EΛ(i)XsXt−EΛ(j)XsXt)(5)

View Source\eqalignno{&{\cal J}(i, j)\cr &\quad ={\cal J}(\Lambda^{(i)}, \Lambda^{(j)})\cr &\quad =\sum_{(s, t)\in E_{i}\cup E_{j}}(\lambda_{st}^{i}-\lambda_{st}^{j})(E_{\Lambda^{(i)}}X_{s}X_{t}-E_{\Lambda^{(j)}}X_{s}X_{t}) &\hbox{(5)} }

Corollary 1:

Let Λ(1) be any Ising model. Let Λ(2) be the model obtained by increasing (decreasing) the parameter corresponding to any edge. The correlation on the edge increases (decreases) in model Λ(2) since J(1,2)≥0. □

Lemma 1:

For all models Λ(i) and Λ(j),D(i,j)≤12J(i,j).

Proof 1:

The lemma follows by noting that

D(i,j)≤J(i,i+j2)+J(j,i+j2),

View Source{\cal D}(i, j)\leq {\cal J}\left(i, {i+j\over 2}\right)+{\cal J}\left(j, {i+j\over 2}\right), and with a little bit of algebra that

J(i,i+j2)+J(j,i+j2)=12J(i,j).□

View Source{\cal J}\left(i, {i+j\over 2}\right)+{\cal J}\left(j, {i+j\over 2}\right)={1\over 2}{\cal J}(i, j).\qquad \square

B. Inference for various values of B

Let Km be a fully connected graph on m+1 vertices. Let Λst be a model with its graph being Kr with the edge (s,t) removed, and every edge parameter being λ. We compute the indirect correlation on the missing edge. In the Ferromagnetic case, the FKG inequality implies that this is the maximum possible indirect correlation due to a model with any graph on these m+1 nodes, and edge parameters being λ.

Lemma 2:

If B≥1, for the model Λst defined above

EΛstXsXt≥1−2(m+1)e3λ/2eB/2+(m+1)e3λ/2

View SourceE_{\Lambda_{{\rm s}t}}X_{s}X_{t}\geq 1-{2(m+1)e^{3\lambda}/2\over e^{B/2}+(m+1)e^{3\lambda/2}}

Proof 2:

A simple computation yields

P(XsXt=1)P(XsXt=−1)=∑ml=0(ml)exp(λ2[(2l−m+1)2−4])∑ml=0(ml)exp(λ2[(2l−m)2])

View Source\eqalignno{&{{\BBP}(X_{s}X_{t}=1)\over {\BBP}(X_{s}X_{t}=-1)}=\cr &\quad{\sum_{l=0}^{m}{ m\choose l }\exp({\lambda\over 2}[(2l-m+1)^{2}-4])\over \sum_{l=0}^{m}{ m\choose l }\exp({\lambda\over 2}[(2l-m)^{2}])}}

We lower bound the above ratio. To do so, we pick the largest term in the denominator. It will be shown that for B≥12logm, the largest term is the one corresponding to i=m and i=0, while for 1≤B≤12logm, the largest term lies in the range i>3m/4 and i<m/4.

Let the largest term be l∗, wolog, ≥m/2. It follows that

P(XsXt=1)P(XsXt=−1)≥(a)(ml∗)exp(λ2[(2l∗−m+1)2−4])(m+1)(ml∗)exp(λ2[(2l∗−m)2])=exp(λ2[4l∗−2m−3])m+1≥exp(λ2[m−3])m+1=exp(B2−32λ)m+1.□

View Source\eqalignno{{{\BBP}(X_{s}X_{t}=1)\over {\BBP}(X_{s}X_{t}=-1)}&\buildrel{(a)}\over {\geq}{{ m\choose l^{\ast} }\exp\left({\lambda\over 2}\left[(2l^{\ast}-m+1)^{2}-4\right]\right)\over (m+1){ m\choose l^{\ast} }\exp\left({\lambda\over 2}\left[(2l^{\ast }-m)^{2}\right]\right)}\cr &={\exp({\lambda\over 2}[4l^{\ast }-2m-3])\over m+1}\cr &\geq{\exp({\lambda\over 2}[m-3])\over m+1}\cr &={\exp({B\over 2}-{3\over 2}\lambda)\over m+1}. \qquad\qquad\qquad\square }

Consider two models in Gd,p,θ with all their parameters being λ:(i)Λ(1) with its graph being Kd with an edge, say (s,t), removed and (ii) Λ(2) with its graph being Kd with a different edge, say (s′,t′), removed. From the FKG inequality and noting that model Λ(2) contains the edge (s,t), observe that

P2(XsXt=1)P2(XsXt=−1)≤P1(XsXt=1)P1(XsXt=−1)e2λ.

View Source{{\BBP}_{2}(X_{s}X_{t}=1)\over {\BBP}_{2}(X_{s}X_{t}=-1)}\leq{{\BBP}_{1}(X_{s}X_{t}=1)\over {\BBP}_{1}(X_{s}X_{t}=-1)}e^{2\lambda}.

Hence

J(1,2)≤2(e2λ−1)P1(XsXt=1)P1(XsXt=−1).(6)

View Source{\cal J}(1,2)\leq{2(e^{2\lambda}-1)\over {{\BBP}_{1}(X_{s}X_{t}=1)\over {\rm P}_{1}(X_{s}X_{t}=-1)}}. \eqno{\hbox{(6)}}

Hence if B>1, it follows that

J(1,2)=λ(EΛ(2)XsXt−EΛ(1)XsXt)+λ(EΛ(1)Xs′Xt′−EΛ(1)Xs′Xt′)≤2Be5λ/2sinh(λ)eB/2.

View Source\eqalignno{{\cal J}(1,2)&=\lambda(E_{\Lambda^{(2)}}X_{s}X_{t}-E_{\Lambda^{(1)}}X_{s}X_{t})+\cr &\qquad\lambda(E_{\Lambda^{(1)}}X_{s^{\prime}}X_{t^{\prime}}-E_{\Lambda^{(1)}}X_{s^{\prime}}X_{t^{\prime}})\cr&\leq{2Be^{5\lambda/2}\sinh(\lambda)\over e^{B/2}}. }

Similarly considering K2k√ and constructing two models from it as above, we obtain two models in Gk,p,θ separated by a distance that is exponentially small in B.

It follows that if B≥1 and any model is to be infered with high probability, the sample complexity should grow exponentially in the parameter B

Lemma 3:

If B≤12, then

EΛstXsXt≤1−2e2B+1

View SourceE_{\Lambda_{{\rm s}t}}X_{s}X_{t}\leq 1-{2\over e^{2B}+1}

Proof 3:

From the definition of B, it follows that

P(XsXt=1)P(XsXt=−1)≤e2B,

View Source{{\BBP}(X_{s}X_{t}=1)\over {\BBP}(X_{s}X_{t}=-1)}\leq e^{2B}, implying the lemma.

We briefly survey some of the results that form the core of the attack on Theorems 1 and 2 in the next two subsections respectively.

B. Necessary number of samples

To estimate the necessary number of samples for decoding the edge structure, we fall back on a version of Fano's inequality. Essentially, Fano's inequality quantifies how much information is gleaned from any single sample, and therefore the number of samples we need before we have sufficient information to infer the edge structure.

1) Fano's inequality

We formalize the above intuition in the following lemma.

Lemma 4:

For any set of w models Λ(1) through Λw and all decoders q:Xn→{Λ(1),…,Λw}, if

n≤w2logw42∑1≤i<j≤wJ(i,j),

View Sourcen\leq{w^{2}\log{w\over 4}\over 2\sum_{1\leq i<j\leq w}{\cal J}(i, j)}, then

max1≤i≤w P(q(Xn)≠Λ(i)|Λ(i))≥12.

View Source{\max_{1\leq i\leq w}\ {\BBP}\left(q(X^{n})\neq\Lambda^{(i)}\vert \Lambda^{(i)}\right)}\geq{1\over 2}.

Proof 4:

Let Λ be a uniform random variable taking values from the set of models Λ(1) through Λw. Let

Pe=1w ∑1≤i≤wP(q(Xn)≠Λ(i)|Λ(i)).

View SourceP_{e}={1\over w}\ \sum\limits_{1\leq i\leq w}{\BBP}\left(q(X^{n})\neq\Lambda^{(i)}\vert \Lambda^{(i)}\right).

From Fano's inequality, H(Λ|Xn)≤1−Pelog(w−1), we obtain

log(w)−1−Pelog(w−1)≤H(Λ)−H(Λ|Xn)=H(Xn)−H(Xn|Λ)=nw∑i=1wD(Λ(i)∥1w∑j=1wΛ(j))≤nw2∑1≤i,j≤wD(i∥j)=nw2 ∑1≤i<j≤w J(i,j).

View Source\eqalignno{\log(w)&-1-P_{e}\log(w-1)\cr &\leq H(\Lambda)-H(\Lambda\vert X^{n})\cr &=H(X^{n})-H(X^{n}\vert \Lambda)\cr &={n\over w}\sum_{i=1}^{w}D(\Lambda^{(i)}\Vert {1\over w}\sum_{j=1}^{w}\Lambda^{(j)})\cr &\leq{n\over w^{2}}\sum_{1\leq i, j\leq w}D(i\Vert j)\cr &={n\over w^{2}}\ \sum_{1\leq {i}<j\leq w}\ {\cal J}(i, j). }

The lemma follows. □

We require a bound on the distance D(i,j) between models Λ(i) and Λ(j) in order to see how far apart the vector of edgewise correlations of the two distributions are.

It is intuitive that the distance D(i,j) between models should be related to the number of edges that exist only in one model. Indeed, Theorem 3 confirms this— D(i,j) is proportional to the matching number of (Ei−Ej)∪(Ej−Ei). Recall that a matching of a graph G is a subgraph H of G such that each vertex in H has degree 1, and that the matching number of G is the largest possible number of edges in any matching.