High-dimensional Ising model selection using ${\ell_1}$-regularized logistic regression

Introduction

Undirected graphical models, also known as Markov random fields, are used in a variety of domains, including statistical physics Ising25 , natural language processing Manning99 , image analysis Woods78 , Hassner80 , Cross83 and spatial statistics Ripley81 , among others. A Markov random field (MRF) is specified by an undirected graph $G=(V,E)$ with vertex set $V=\{1,2,\ldots,p\}$ and edge set $E\subset V\times V$. The structure of this graph encodes certain conditional independence assumptions among subsets of the $p$-dimensional discrete random variable $X=(X_{1},X_{2},\ldots,X_{p})$ where variable $X_{i}$ is associated with vertex $i\in V$. One important problem for such models is to estimate the underlying graph from $n$ independent and identically distributed samples $\{x^{(1)},x^{(2)},\ldots,x^{(n)}\}$ drawn from the distribution specified by some Markov random field. As a concrete illustration, for binary random variables, each vector-valued sample $x^{(i)}\in\{0,1\}^{p}$ might correspond to the votes of a set of $p$ politicians on a particular bill, and estimating the graph structure amounts to detecting statistical dependencies in these voting patterns (see Banerjee, Ghaoui and d’Asprémont BanGhaAsp08 for further discussion of this example).

Due to both its importance and difficulty, the problem of structure learning for discrete graphical models has attracted considerable attention. The absence of an edge in a graphical model encodes a conditional independence assumption. Constraint-based approaches spirtes00 estimate these conditional independencies from the data using hypothesis testing and then determine a graph that most closely represents those independencies. Each graph represents a model class of graphical models; learning a graph then is a model class selection problem. Score-based approaches combine a metric for the complexity of the graph with a measure of the goodness of fit of the graph to the data; for instance, log-likelihood of the maximum likelihood parameters given the graph, to obtain a score for each graph. The score is used together with a search procedure that generates candidate graph structures to be scored. The number of graph structures grows super-exponentially, however, and Chickering chickering95 shows that this problem is in general NP-hard.

A complication for undirected graphical models involving discrete random variables is that typical score metrics involve the partition function or cumulant function associated with the Markov random field. For general undirected MRFs, calculation of this partition function is computationally intractable Welsh93 . The space of candidate structures in scoring based approaches is thus typically restricted to either directed graphical models dasgupta99 or to simple sub-classes of undirected graphical models such as those based on trees chowliu68 and hypertrees srebro03 . Abbeel, Koller and Ng AbbKolNg06 propose a method for learning factor graphs based on local conditional entropies and thresholding and analyze its behavior in terms of Kullback–Leibler divergence between the fitted and true models. They obtain a sample complexity that grows logarithmically in the number of vertices $p$, but the computational complexity grows at least as quickly as ${\mathcal{O}}(p^{d+1})$ where $d$ is the maximum neighborhood size in the graphical model. This order of complexity arises from the fact that for each node, there are ${p\choose d}={\mathcal{O}}(p^{d})$ possible neighborhoods of size $d$ for a graph with $p$ vertices. Csiszár and Talata Csiszar06 show consistency of a method that uses pseudo-likelihood and a modification of the BIC criterion, but this also involves a prohibitively expensive search.

Portions of this work were initially reported in a conference publication WaiRavLaf06 , with the weaker result that $n=\Omega(d^{6}\log d+d^{5}\log p)$ samples suffice for consistent Ising model selection. Since the appearance of that paper, other researchers have also studied the problem of model selection in discrete Markov random fields. For the special case of bounded degree models, Bresler, Mossel and Sly Bresler08 describe a simple search-based method, and prove under relatively mild assumptions that it can recover the graph structure with $\Theta(\log p)$ samples. However, in the absence of additional restrictions, the computational complexity of the method is ${\mathcal{O}}(p^{d+1})$. In other work, Santhanam and Wainwright SanWai08 analyze the information-theoretic limits of graphical model selection, providing both upper and lower bounds on various model selection procedures, but these methods also have prohibitive computational costs.

The remainder of this paper is organized as follows. We begin in Section 2 with background on discrete graphical models, the model selection problem and logistic regression. In Section 3, we state our main result, develop some of its consequences and provide a high-level outline of the proof. Section 4 is devoted to proving a result under stronger assumptions on the sample Fisher information matrix whereas Section 5 provides concentration results linking the population matrices to the sample versions. In Section 6, we provide some experimental results that illustrate the practical performance of our method and the close agreement between theory and practice. Section 7 discusses an extension to more general Markov random fields, and we conclude in Section 8.

Background and problem formulation

We begin by providing some background on Markov random fields, defining the problem of graphical model selection and describing our method based on neighborhood logistic regression.

The partition function $Z({{\theta^{*}}})$ ensures that the distribution sums to one. This model is used in many applications of spatial statistics such as modeling the behavior of gases or magnets in statistical physics Ising25 , building statistical models in computer vision Geman84 and social network analysis.

2 Graphical model selection

Here the sign function takes value $+1$ if $\theta^{*}_{st}>0$, value $-1$ if $\theta^{*}_{st}<0$ and , otherwise. Note that the weaker graphical model selection problem amounts to recovering the vector $|E^{*}|$ of absolute values.

The classical notion of statistical consistency applies to the limiting behavior of an estimation procedure as the sample size $n$ goes to infinity with the model size $p$ itself remaining fixed. In many contemporary applications of graphical models—among them gene microarray data and social network analysis—the model dimension $p$ is comparable to or larger than the sample size $n$, so that the relevance of such “fixed $p$” asymptotics is limited. With this motivation, our analysis in this paper is of the high-dimensional nature, in which both the model dimension and the sample size are allowed to increase, and we study the scalings under which consistent model selection is achievable.

More precisely, we consider sequences of graphical model selection problems, indexed by the sample size $n$, number of vertices $p$ and maximum node degree $d$. We assume that the sample size $n$ goes to infinity, and both the problem dimension $p=p(n)$ and $d=d(n)$ may also scale as a function of $n$. The setting of fixed $p$ or $d$ is covered as a special case. Let $\widehat{E}_{n}$ be an estimator of the signed edge pattern $E^{*}$ based on the $n$ samples. Our goal is to establish sufficient conditions on the scaling of the triple $(n,p,d)$ such that our proposed estimator is consistent in the sense that

We sometimes call this property sparsistency, as a shorthand for consistency of the sparsity pattern of the parameter $\theta^{*}$.

3 Neighborhood-based logistic regression

Recovering the signed edge vector $E^{*}$ of an undirected graph $G$ is equivalent to recovering, for each vertex $r\in V$, its neighborhood set $\mathcal{N}(r):=\{t\in V\mid(r,t)\in E\}$ along with the correct signs $\operatorname{sign}(\theta^{*}_{rt})$ for all $t\in\mathcal{N}(r)$. To capture both the neighborhood structure and sign pattern, we define the product set of “signed vertices” as $\{-1,1\}\times V$. We use the shorthand “$\iota r$” for elements $(\iota,r)\in\{-1,1\}\times V$. We then define the signed neighborhood set as

Here the sign function has an unambiguous definition, since $\theta^{*}_{rt}\neq 0$ for all $t\in\mathcal{N}(r)$. Observe that this signed neighborhood set $\mathcal{N}_{\pm}(r)$ can be recovered from the sign-sparsity pattern of the $(p-1)$-dimensional subvector of parameters

associated with vertex $r$. In order to estimate this vector $\theta^{*}_{\setminus r}$, we consider the structure of the conditional distribution of $X_{r}$ given the other variables $X_{\setminus r}=\{X_{t}\mid t\in V\setminus\{r\}\}$. A simple calculation shows that under the model (1), this conditional distribution takes the form

Thus the variable $X_{r}$ can be viewed as the response variable in a logistic regression in which all of the other variables $X_{\setminus r}$ play the role of the covariates.

is the rescaled negative log likelihood (the rescaling factor $1/n$ in this definition is for later theoretical convenience) and $\lambda_{(n,p,d)}>0$ is a regularization parameter, to be specified by the user. For notational convenience, we will also use $\lambda_{n}$ as notation for this regularization parameter suppressing the potential dependence on $p$ and $d$.

Following some algebraic manipulation, the regularized negative log likelihood can be written as

Accordingly, let $\widehat{\theta}^{n}_{\setminus r}$ be an element of the minimizing set of problem (7). Although $\widehat{\theta}^{n}_{\setminus r}$ need not be unique in general since the problem (7) need not be strictly convex, our analysis shows that in the regime of interest, this minimizer $\widehat{\theta}^{n}_{\setminus r}$ is indeed unique. We use $\widehat{\theta}^{n}_{\setminus r}$ to estimate the signed neighborhood $\mathcal{N}_{\pm}(r)$ according to

We say that the full graph $G$ is estimated consistently, written as the event $\{\widehat{E}_{n}=E^{*}\}$, if every signed neighborhood is recovered—that is, $\widehat{\mathcal{N}}_{\pm}(r)=\mathcal{N}_{\pm}(r)$ for all $r\in V$.

Method and theoretical guarantees

Our main result concerns conditions on the sample size $n$ relative to the parameters of the graphical model—more specifically, the number of nodes $p$ and maximum node degree $d$—that ensure that the collection of signed neighborhood estimates (9), one for each node $r$ of the graph, agree with the true neighborhoods so that the full graph is estimated consistently. In this section, we begin by stating the assumptions that underlie our analysis, and then give a precise statement of the main result. We then provide a high-level overview of the key steps involved in its proof, deferring details to later sections. Our analysis proceeds by first establishing sufficient conditions for correct signed neighborhood recovery—that is, $\{\widehat{\mathcal{N}}_{\pm}(r)=\mathcal{N}_{\pm}(r)\}$—for some fixed node $r\in V$. By showing that this neighborhood consistency is achieved at sufficiently fast rates, we can then use a union bound over all $p$ nodes of the graph to conclude that consistent graph selection is also achieved.

For future reference, this is given as the explicit expression

In the following we write simply $Q^{*}$ for the matrix $Q^{*}_{r}$ where the reference node $r$ should be understood implicitly. Moreover, we use $S:=\{(r,t)\mid t\in\mathcal{N}(r)\}$ to denote the subset of indices associated with edges of $r$, and ${{S^{c}}}$ to denote its complement. We use $Q^{*}_{SS}$ to denote the $d\times d$ sub-matrix of $Q^{*}$ indexed by $S$. With this notation, we state our assumptions:

The subset of the Fisher information matrix corresponding to the relevant covariates has bounded eigenvalues; that is, there exists a constant $C_{\min}>0$ such that

(A2) Incoherence condition

Our next assumption captures the intuition that the large number of irrelevant covariates (i.e., nonneighbors of node $r$) cannot exert an overly strong effect on the subset of relevant covariates (i.e., neighbors of node $r$). To formalize this intuition, we require the existence of an $\alpha\in(0,1]$ such that

2 Statement of main result

We are now ready to state our main result on the performance of neighborhood logistic regression for graphical model selection. Naturally, the limits of model selection are determined by the minimum value over the parameters $\theta^{*}_{rt}$ for pairs $(r,t)$ included in the edge set of the true graph. Accordingly, we define the parameter

With this definition, we have the following:

Consider an Ising graphical model with parameter vector $\theta^{*}$ and associated edge set $E^{*}$ such that conditions (A1) and (A2) are satisfied by the population Fisher information matrix $Q^{*}$, and let $\mathfrak{X}_{1}^{n}$ be a set of $n$ i.i.d. samples from the model specified by $\theta^{*}$. Suppose that the regularization parameter $\lambda_{n}$ is selected to satisfy

Then there exist positive constants $L$ and $K$, independent of $(n,p,d)$, such that if

then the following properties hold with probability at least $1-2\exp(-K\lambda_{n}^{2}n)$.

For each $r\in V$, the estimated signed neighborhood $\widehat{\mathcal{N}}_{\pm}(r)$ correctly excludes all edges not in the true neighborhood. Moreover, it correctly includes all edges $(r,t)$ for which $|\theta_{rt}^{*}|\geq\frac{10}{C_{\min}}\sqrt{d}\lambda_{n}$.

The theorem not only specifies sufficient conditions but also the probability with which the method recovers the true signed edge-set. This probability decays exponentially as a function of $\lambda_{n}^{2}n$ which leads naturally to the following corollary on model selection consistency of the method for a sequence of Ising models specified by $(n,p(n),d(n))$.

Consider a sequence of Ising models with graph edge sets $\{E^{*}_{p(n)}\}$ and parameters $\{\theta^{*}_{(n,p,d)}\}$; each of which satisfies conditions (A1) and (A2). For each $n$, let $\mathfrak{X}_{1}^{n}$ be a set of $n$ i.i.d. samples from the model specified by $\theta^{*}_{(n,p,d)}$, and suppose that $(n,p(n),d(n))$ satisfies the scaling condition (17) of Theorem 1. Suppose further that the sequence $\{\lambda_{n}\}$ of regularization parameters satisfies condition (16) and

and the minimum parameter weights satisfy

(a) It is worth noting that the scaling condition (17) on $(n,p,d)$ allows for graphs and sample sizes in the “large $p$, small $n$” regime (meaning $p\gg n$), as long as the degrees are bounded, or grow at a sufficiently slow rate. In particular, one set of sufficient conditions are the scalings

for some constants $c_{1},c_{2}>0$. Under these scalings, note that we have $d^{3}\log(p)=\mathcal{O}(n^{3c_{1}+c_{2}})=o(n)$, so that condition (17) holds.

A bit more generally, note that in the regime $p\gg n$, the growth condition (17) requires that that $d=o(p)$. However, in many practical applications of graphical models (e.g., image analysis, social networks), one is interested in node degrees $d$ that remain bounded or grow sub-linearly in the graph size so that this condition is not unreasonable.

(b) Loosely stated, the theorem requires that the edge weights are not too close to zero (in absolute value) for the method to estimate the true graph. In particular, conditions (16) and (19) imply that the minimum edge weight $\theta^{*}_{\min}$ is required to scale as

Note that in the classical fixed $(p,d)$ case, this reduces to the familiar scaling requirement of $\theta^{*}_{\min}=\Omega(n^{-1/2})$.

(c) In the high-dimensional setting (for $p\rightarrow+\infty)$, a choice of the regularization parameter satisfying both conditions (16) and (18) is, for example,

for which the probability of incorrect model selection decays at rate${\mathcal{O}}(\exp(-K^{\prime}\log p))$ for some constant $K^{\prime}>0$. In the classical setting (fixed $p$), this choice can be modified to $\lambda_{n}=\frac{16(2-\alpha)}{\alpha}\sqrt{\frac{\log(pn)}{n}}$.

The analysis required to prove Theorem 1 can be divided naturally into two parts. First, in Section 4, we prove a result (stated as Proposition 1) for “fixed design” matrices. More precisely, we show that if the dependence condition (A1) and the mutual incoherence condition (A2) hold for the sample Fisher information matrix

then the growth condition (17) and choice of $\lambda_{n}$ from Theorem 1 are sufficient to ensure that the graph is recovered with high probability.

The second part of the analysis, provided in Section 5, is devoted to showing that under the specified growth condition (17), imposing incoherence and dependence assumptions on the population version of the Fisher information $Q^{*}$ guarantees (with high probability) that analogous conditions hold for the sample quantities $Q^{n}$. On one hand, it follows immediately from the law of large numbers that the empirical Fisher information $Q^{n}_{AA}$ converges to the population version $Q^{*}_{AA}$ for any fixed subset $A$. However, in the current setting, the added delicacy is that we are required to control this convergence over subsets of increasing size. Our proof therefore requires some large-deviation analysis for random matrices with dependent elements so as to provide exponential control on the rates of convergence.

3 Primal-dual witness for graph recovery

For the convex program (7), the zero sub-gradient optimality conditions Rockafellar take the form

Based on this lemma, we construct a primal-dual witness $(\widehat{\theta},\widehat{z})$ with the following steps.

First, we set $\widehat{\theta}_{S}$ as the minimizer of the partial penalized likelihood

and set $\widehat{z}_{S}=\operatorname{sign}(\widehat{\theta}_{S})$.

Second, we set $\widehat{\theta}_{{S^{c}}}=0$ so that condition (23b) holds.

In the third step, we obtain $\widehat{z}_{{{S^{c}}}}$ from (21) by substituting in the values of $\widehat{\theta}$ and $\widehat{z}_{S}$. Thus our construction satisfies conditions (23b) and (21).

The final and most challenging step consists of showing that the stated scalings of $(n,p,d)$ imply that, with high-probability, the remaining conditions (23a) and (22) are satisfied.

Our analysis in step (d) guarantees that $\|\widehat{z}_{S^{c}}\|_{\infty}<1$ with high probability. Moreover, under the conditions of Theorem 1, we prove that the sub-matrix of the sample Fisher information matrix is strictly positive definite with high probability so that by Lemma 1, the primal solution $\widehat{\theta}$ is guaranteed to be unique.

Analysis under sample Fisher matrix assumptions

As in the statement of Theorem 1, the quantities $L$ and $K$ refer to constants independent of $(n,p,d)$. With this notation, we have the following:

If the event $\mathcal{M}(\mathfrak{X}_{1}^{n})$ holds, the sample size satisfies $n>Ld^{2}\log(p)$, and the regularization parameter is chosen such that $\lambda_{n}\geq\frac{16(2-\alpha)}{\alpha}\sqrt{\frac{\log p}{n}}$, then with probability at least $1-2\exp(-K\lambda_{n}^{2}n)\rightarrow 1$, the following properties hold.

(b) For each $r\in V$, the estimated signed neighborhood vector $\widehat{\mathcal{N}}_{\pm}(r)$ correctly excludes all edges not in the true neighborhood. Moreover, it correctly includes all edges with $|\theta_{rt}|\geq\frac{10}{C_{\min}}\sqrt{d}\lambda_{n}$.

Loosely stated, this result guarantees that if the sample Fisher information matrix is “good,” then the conditional probability of successful graph recovery converges to zero at the specified rate. The remainder of this section is devoted to the proof of Proposition 1.

We begin with statements of some key technical lemmas that are central to our main argument with their proofs deferred to Appendix B. The central object is the following expansion obtained by re-writing the zero-subgradient condition as

with $\accentset{\vspace*{-2pt}-}{\theta}^{(j)}$ a parameter vector on the line between ${{\theta^{*}}}$ and $\widehat{\theta}$, and with $[\cdot]_{j}^{T}$ denoting the $j$th row of the matrix. The following lemma addresses the behavior of the term $W^{n}$ in this expansion:

For the specified mutual incoherence parameter $\alpha\in(0,1]$, we have

which converges to zero at rate $\exp(-c\lambda_{n}^{2}n)$ as long as $\lambda_{n}\geq\frac{16(2-\alpha)}{\alpha}\sqrt{\frac{\log p}{n}}$.

See Appendix B.1 for the proof of this claim.

If $\lambda_{n}d\leq\frac{C_{\min}^{2}}{10D_{\max}}$ and$\|W^{n}\|_{\infty}\leq\lambda_{n}/4$, then

See Appendix B.2 for the proof of this claim.

Our final technical lemma provides control on the remainder term (28).

If $\lambda_{n}d\leq\frac{C_{\min}^{2}}{100D_{\max}}\frac{\alpha}{2-\alpha}$ and $\|W^{n}\|_{\infty}\leq\lambda_{n}/4$, then

See Appendix B.3 for the proof of this claim.

2 Proof of Proposition 1

Using these lemmas, the proof of Proposition 1 is straightforward. Consider the choice of the regularization parameter, $\lambda_{n}=16\frac{2-\alpha}{\alpha}\sqrt{\frac{\log p}{n}}$. This choice satisfies the condition of Lemma 2, so that we may conclude that with probability greater than $1-2\exp(-c\lambda_{n}^{2}n)\rightarrow 1$, we have

using the fact that $\alpha\leq 1$. The remaining two conditions that we need to apply the technical lemmas concern upper bounds on the quantity $\lambda_{n}d$. In particular, for a sample size satisfying $n>\frac{100^{2}D_{\max}^{2}}{C_{\min}^{4}}\frac{(2-\alpha)^{4}}{\alpha^{4}}d^{2}\log p$, we have

so that the conditions of both Lemmas 3 and 4 are satisfied.

Since the matrix $Q^{n}_{SS}$ is invertible by assumption, the conditions (4.2) can be re-written as

We now demonstrate that for the dual sub-vector $\widehat{z}_{S^{c}}$ defined by (33), we have $\|\widehat{z}_{S^{c}}\|_{\infty}<1$. Using the triangle inequality and the mutual incoherence bound (14), we have that

Correct sign recovery

We next show that our primal sub-vector $\widehat{\theta}_{S}$ defined by (24) satisfies sign consistency, meaning that $\operatorname{sgn}(\widehat{\theta}_{S})=\operatorname{sgn}(\theta^{*}_{S})$. In order to do so, it suffices to show that

recalling the notation $\theta^{*}_{\min}:=\min_{(r,t)\in E}|\theta^{*}_{rt}|$. From Lemma 3, we have $\|\theta_{S}-\theta^{*}_{S}\|_{2}\leq\frac{5}{C_{\min}}\sqrt{d}\lambda_{n}$ so that

which is less than one as long as $\theta^{*}_{\min}\geq\frac{10}{C_{\min}}\sqrt{d}\lambda_{n}$.

Uniform convergence of sample information matrices

In this section we complete the proof of Theorem 1 by showing that if the dependency (A1) and incoherence (A2) assumptions are imposed on the population Fisher information matrix then under the specified scaling of $(n,p,d)$, analogous bounds hold for the sample Fisher information matrices with probability converging to one. These results are not immediate consequences of classical random matrix theory (e.g., DavSza01 ) since the elements of $Q^{n}$ are highly dependent. Recall the definitions

The following result is the analog for the incoherence assumption (A2) showing that the scaling of $(n,p,d)$ given in Theorem 1 guarantees that population incoherence implies sample incoherence.

If the population covariance satisfies a mutual incoherence condition (14) with parameter $\alpha\in(0,1]$ as in assumption (A2), then the sample matrix satisfies an analogous version, with high probability in the sense that

Proofs of these two lemmas are provided in the following sections. Before proceeding, we take note of a simple bound to be used repeatedly throughout our arguments. By definition of the matrices $Q^{n}(\theta)$ and $Q(\theta)$ [see (20) and (11)], the $(j,k)$th element of the difference matrix $Q^{n}(\theta)-Q(\theta)$ can be written as an i.i.d. sum of the form $Z_{jk}=\frac{1}{n}\sum_{i=1}^{n}Z^{(i)}_{jk}$ where each $Z^{(i)}_{jk}$ is zero-mean and bounded (in particular, $|Z^{(i)}_{jk}|\leq 4$). By the Azuma–Hoeffding bound Hoeffding63 , for any indices $j,k=1,\ldots,d$ and for any $\varepsilon>0$, we have

So as to simplify notation, throughout this section, we use $K$ to denote a universal positive constant, independent of $(n,p,d)$. Note that the precise value and meaning of $K$ may differ from line to line.

By the Courant–Fischer variational representation Horn85 , we have

Hence it suffices to obtain a bound on the spectral norm $|\!|\!|Q_{SS}-Q^{n}_{SS}|\!|\!|_{{2}}$. Observe that

Setting $\varepsilon^{2}=\delta^{2}/d^{2}$ in (39) and applying the union bound over the $d^{2}$ index pairs $(j,k)$ then yields

which obeys the same upper bound (40) by following the analogous argument.

2 Proof of Lemma 6

We begin by decomposing the sample matrix as the sum $Q^{n}_{{S^{c}}S}(Q^{n}_{SS})^{-1}=T_{1}+T_{2}+T_{3}+T_{4}$ where we define

The fourth term is easily controlled; indeed, we have

by the incoherence assumption (A2). If we can show that $|\!|\!|T_{i}|\!|\!|_{{\infty}}\leq\frac{\alpha}{6}$ for the remaining indices $i=1,2,3$, then by our four term decomposition and the triangle inequality, the sample version satisfies the bound (38), as claimed. We deal with these remaining terms using the following lemmas:

For any $\delta>0$ and constants $K,K^{\prime}$, the following bounds hold:

See Appendix C for the proof of these claims.

Turning to the first term, we re-factorize it as

and then bound it (using the sub-multiplicative property $|\!|\!|AB|\!|\!|_{{\infty}}\leq\break|\!|\!|A|\!|\!|_{{\infty}}|\!|\!|B|\!|\!|_{{\infty}}$) as follows:

where we have used the incoherence assumption (A2). Using the bound (37b) from Lemma 5 with $\delta=C_{\min}/2$, we have $|\!|\!|(Q^{n}_{SS})^{-1}|\!|\!|_{{2}}=[\Lambda_{\min}(Q^{n}_{SS})]^{-1}\leq\frac{2}{C_{\min}}$ with probability greater than $1-\exp(-Kn/d^{2}+2\log(d))$. Next, applying the bound (7) with $\delta=c/\sqrt{d}$, we conclude that with probability greater than $1-2\exp(-Knc^{2}/d^{3}+\log(d))$, we have

By choosing the constant $c>0$ sufficiently small, we are guaranteed that

Control of second term

We then apply bound (7) with $\delta=\frac{\alpha}{3}\frac{C_{\min}}{\sqrt{d}}$ to conclude that

Control of third term

Finally, in order to bound the third term $T_{3}$, we apply the bounds (7) and (7), both with $\delta=\sqrt{\alpha/3}$, and use the fact that $\log(d)\leq\log(p-d)$ to conclude that

Putting together all of the pieces, we conclude that

Experimental results

Figure 2 shows results for the $4$-nearest-neighbor grid model, illustrated in Figure 1(a) for three different graph sizes $p\in\{64,100,225\}$ with mixed couplings [panel (a)] and attractive couplings [panel (b)]. Each curve corresponds to a given problem size, and corresponds to the success probability versus the control parameter $\beta$. Each point corresponds to the average of $N=200$ trials. Notice how, despite the very different regimes of $(n,p)$ that underlie each curve, the different curves all line up with one another quite well. This fact shows that for a fixed degree graph (in this case $\deg=4$), the ratio $n/\log(p)$ controls the success/failure of our model selection procedure which is consistent with the prediction of Theorem 1. Figure 3 shows analogous results for the $8$-nearest-neighbor lattice model ($d=8$), for the same range of problem size $p\in\{64,100,225\}$ and for both mixed and attractive couplings. Notice how once again the curves for different problem sizes are all well aligned which is consistent with the prediction of Theorem 1.

For our next set of experiments, we investigate the performance of our method for a class of graphs with unbounded maximum degree $d$. In particular, we construct star-shaped graphs with $p$ vertices by designating one node as the hub and connecting it to $d<(p-1)$ of its neighbors. For linear sparsity, we choose $d=\lceil 0.1p\rceil$, whereas for logarithmic sparsity we choose $d=\lceil\log(p)\rceil$. We again study a triple of graph sizes $p\in\{64,100,225\}$, and Figure 4 shows the resulting curves of success probability versus control parameter $\beta=n/[10d\log(p)]$. Panels (a) and (b) correspond, respectively, to the cases of logarithmic and linear degrees. As with the bounded degree models in Figure 2 and 3, these curves align with one another showing a transition from failure to success with probability one.

For comparative purposes, we also illustrate the performance of the PC algorithm of Spirtes, Glymour and Scheines spirtes00 as well as the maximum weight tree method of Chow and Liu chowliu68 . Since the star graph is a tree (cycle-free), both of these methods are applicable in this case. The PC algorithm is targeted to learning (equivalence classes of) directed acyclic graphs, and consists of two stages. In the first stage it starts from a completely connected undirected graph, and iteratively removes edges based on conditional independence tests so that at the end of this stage it is left with an undirected graph which is called a skeleton. In the second stage, it partially directs some of the edges in the skeleton so as to obtain a completed partially directed acyclic graph which corresponds to an equivalence class of directed acyclic graphs. As pointed out by Kalisch and Bühlmann kalisch07 , for high-dimensional problems, the output of the first stage, which is the undirected skeleton graph, could provide a useful characterization of the dependencies in the data. Following this suggestion, we use the skeleton graph determined by the first stage of the PC algorithm as an estimate of the graph structure. We use the pcalg R-package kalisch07 as an implementation of the PC algorithm which uses partial correlations to test conditional independencies.

The Chow–Liu algorithm chowliu68 is a method for exact maximum likelihood structure selection which is applicable to the case of trees. More specifically, it chooses, from among all trees with a specified number of edges, the tree that minimizes the Kullback–Leibler divergence to the empirical distribution defined by the samples. From an implementational point of view, it starts with a completely connected weighted graph with edge weights equal to the empirical mutual information between the incident node variables of the edge and then computes its maximum weight spanning tree. Since our underlying model is a star-shaped graph with fewer than $(p-1)$ edges, a spanning tree would necessarily include false positives. We thus estimate the maximum weight forest with $d$ edges instead where we supplied the number of edges $d$ in the true graph to the algorithm.

Extensions to general discrete Markov random fields

Since $\mathcal{X}$ is discrete, each potential function $\phi_{st}$ can be parameterized as linear combinations of $\{0,1\}$-valued indicator functions. In particular, for each $s\in V$ and $j\in\{1,\ldots,{m}-1\}$, we define

Any set of potential functions can then be written as

With this set-up, we now describe a graph selection procedure that is the natural generalization of our procedure for the Ising model. As before we focus on recovering for each vertex $r\in V$ its neighborhood set and then combine the neighborhood sets across vertices to form the graph estimate.

In order to estimate this column support, we consider the conditional distribution of $X_{r}$ given the other variables $X_{\setminus\{r\}}=\{X_{t}|t\in V\setminus\{r\}\}$. For a binary model, this distribution is of the logistic form while for a general pairwise MRF, it takes the form

Thus, $X_{r}$ can be viewed as the response variable in a multiclass logistic regression in which the indicator functions associated with the other variables,

More specifically, let $\mathfrak{X}_{1}^{n}=\{x^{(1)},\ldots,x^{(n)}\}$ denote an i.i.d. set of $n$ samples, drawn from the discrete MRF (46). In order to estimate the neighborhood of node $r$, we solve the following convex program:

Conclusion

As discussed in Section 7, although the current analysis is applied to binary Markov random fields, the methods of this paper can be extended to general discrete graphical models with a higher number of states using a multinomial likelihood and some form of block regularization. It should also be possible and would be interesting to obtain high-dimensional rates in this setting. A final interesting direction for future work is the case of samples drawn in a non-i.i.d. manner from some unknown Markov random field; we suspect that similar results would hold for weakly dependent sampling schemes.

Appendix A Proof of uniqueness lemma

In this appendix, we prove Lemma 1. By Lagrangian duality, the penalized problem (7) can be written as an equivalent constrained optimization problem over the ball $\|\theta\|_{1}\leq C(\lambda_{n})$, for some constant $C(\lambda_{n})<+\infty$. Since the Lagrange multiplier associated with this constraint—namely, $\lambda_{n}$—is strictly positive, the constraint is active at any optimal solution so that $\|\theta\|_{1}$ is constant across all optimal solutions.

where the weights $\alpha_{v}$ are nonnegative and sum to one. In fact, these weights correspond to an optimal vector of Lagrange multipliers for an alternative formulation of the problem in which $\alpha_{v}$ is the Lagrange multiplier for the constraint $\langle v,\theta\rangle\leq C(\lambda_{n})$. From standard Lagrangian theory Bertsekasnonlin , it follows that any other optimal primal solution $\widetilde{\theta}$ must minimize the associated Lagrangian—or equivalently, satisfy (21)—and moreover must satisfy the complementary slackness conditions $\alpha_{v}\{\langle v,\widetilde{\theta}\rangle-C\}=0$ for all sign vectors $v$. But these conditions imply that $\langle\widehat{z},\widetilde{\theta}\rangle=C=\|\widetilde{\theta}\|_{1}$ which cannot occur if $\widetilde{\theta}_{j}\neq 0$ for some index $j$ for which $|\widehat{z}_{j}|<1$. We thus conclude that $\widetilde{\theta}_{{{S^{c}}}}=0$ for all optimal primal solutions.

Finally, given that all optimal solutions satisfy $\theta_{{{S^{c}}}}=0$, we may consider the restricted optimization problem subject to this set of constraints. If the principal submatrix of the Hessian is positive definite, then this sub-problem is strictly convex so that the optimal solution must be unique.

Appendix B Proofs for technical lemmas

In this section, we provide proofs of Lemmas 2, 3 and 4, previously stated in Section 4.

Note that any entry of $W^{n}$ has the form $W^{n}_{u}=\frac{1}{n}\sum_{i=1}^{n}Z^{(i)}_{u}$ where for $i=1,2,\ldots,n$, the variables

Finally, applying a union bound over the indices $u$ of $W^{n}$ yields

B.2 Proof of Lemma 3

contradicting the assumed strict positivity of $G$ on the boundary.

for some $\alpha\in$. For the first term, we have the bound

since $\|W_{S}\|_{\infty}\leq\frac{\lambda_{n}}{4}$ by assumption.

Applying the triangle inequality to the last term in the expansion (51) yields

Since $\|u_{S}\|_{1}\leq\sqrt{d}\|u_{S}\|_{2}$, we have

Finally, turning to the middle Hessian term, we have

By a Taylor series expansion of $\eta(x^{(i)};\cdot)$, we have

Now note that $|\eta^{\prime}(x^{(i)};\theta^{*}_{S}+\alpha u_{S})|\leq 1$, and

by assumption. Combining these pieces, we obtain

assuming that $\lambda_{n}\leq\frac{C_{\min}}{2MD_{\max}d}$. We verify this condition momentarily, after we have specified the constant $M$.

Finally, combining the bounds (52), (53), and (54) in the expression (51), we conclude that

This expression is strictly positive for $M=5/C_{\min}$. Consequently, as long as

as assumed in the statement of the lemma, we are guaranteed that

B.3 Proof of Lemma 4

We first show that the remainder term $R^{n}$ satisfies the bound $\|R^{n}\|_{\infty}\leq D_{\max}\|\widehat{\theta}_{S}-\theta^{*}_{S}\|_{2}^{2}$. Then the result of Lemma 3—namely, that $\|\widehat{\theta}_{S}-\theta^{*}_{S}\|_{2}\leq\frac{5}{C_{\min}}\lambda_{n}\sqrt{d}$—can be used to conclude that

Focusing on element $R^{n}_{j}$ for some index $j\in\{1,\ldots,p\}$, we have

for some point $\accentset{\vspace*{-2pt}-}{\theta}^{(j)}=t_{j}\widehat{\theta}+(1-t_{j}){{\theta^{*}}}$. Setting $g(t)=\frac{4\exp(2t)}{[1+\exp(2t)]^{2}}$, note that $\eta(x;\theta)=g(x_{r}\sum_{t\in V\setminus r}\theta_{rt}x_{t})$. By the chain rule and another application of the mean value theorem, we then have

where $\accentset{\vspace*{-1pt}=}{\theta}^{(j)}$ is another point on the line joining $\widehat{\theta}$ and ${{\theta^{*}}}$. Setting $a_{i}:=\{g^{\prime}(\accentset{\vspace*{-1pt}=}{\theta}^{(j)T}\times x^{(i)})x^{(i)}_{j}\}$ and $b_{i}:=\{[\accentset{\vspace*{-2pt}-}{\theta}^{(j)}-\theta^{*}]^{T}x^{(i)}(x^{(i)})^{T}[\widehat{\theta}-\theta^{*}]\}$, we have

A calculation shows that $\|a\|_{\infty}\leq 1$, and

where the second line uses the fact that $\widehat{\theta}_{S^{c}}=\theta^{*}_{S^{c}}=0$. This concludes the proof.

Appendix C Proof of Lemma 7

where the final inequality uses a union bound, and the fact that $|{S^{c}}|\leq p-d$. Via another union bound over the row elements

from which the claim (7) follows by setting $\varepsilon=\delta/d$ in the Hoeffding bound (39). The proof of bound (7) is analogous with the pre-factor $(p-d)$ replaced by $d$.

From the proof of Lemma 5, in particular equation (40), we have

for a constant $B$. Moreover, from (40), we have

References