High-Dimensional Graphical Model Selection Using $\ell_1$-Regularized Logistic Regression

Introduction

Undirected graphical models, also known as Markov random fields (MRFs), are used in a variety of domains, including artificial intelligence, natural language processing, image analysis, statistical physics, and spatial statistics, 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$. A fundamental problem is the graphical model selection problem: given a set of $n$ samples $\{x^{(1)},x^{(2)},\ldots,x^{(n)}\}$ from a Markov random field, estimate the structure of the underlying graph. The sample complexity of such an estimator is the minimal number of samples $n$, as a function of the graph size $p$ and possibly other parameters such as the maximum node degree $d$, required for the probability of correct identification of the graph to converge to one. Another important property of any model selection procedure is its computational complexity.

Due to both its importance and difficulty, structure learning in random fields has attracted considerable attention. The absence of an edge in a graphical model encodes a conditional independence assumption. Constraint-based approaches (Spirtes et al. 2000) 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 goodness of fit measure 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 then 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 1995 shows that this problem is in general NP-hard.

A complication for undirected graphical models is that typical score metrics involve the normalization constant (also called the partition function) associated with the Markov random field, which is intractable (#P) to compute for general undirected models. The space of candidate structures in scoring based approaches is thus typically restricted to either directed models—Bayesian networks—or simple undirected graph classes such as trees (Chow and Liu 1968), polytrees (Dasgupta 1999) and hypertrees (Srebro 2003). Abbeel et al. 2006 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 graph size $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 $\binom{p}{d}={\mathcal{O}}(p^{d})$ possible neighborhoods of size $d$ for a graph with $p$ vertices. Csiszár and Talata 2006 show consistency of a method that uses pseudo-likelihood and a modification of the BIC criterion, but this also involves a prohibitively expensive search.

In work subsequent to the initial conference version of this work (Wainwright et al. 2007), other researchers have also studied the problem of model selection in discrete Markov random fields. For the special case of bounded degree models, Bresler et al. 2008 describe a simple search-based method, and prove under relatively mild assumptions that it can recover the graph structure with $\Theta(d\log p)$ samples. However, in the absence of additional restrictions, the computational complexity of the method is ${\mathcal{O}}(p^{d+1})$. Santhanam and Wainwright 2008 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 cost.

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 itself, whereas Section 5 provides concentration results linking the population matrices to the sample versions. In Section 6, we provide some experimental results to illustrate the practical performance of our method, and the close agreement between theory and practice, and we conclude in Section 7.

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. The Ising model has proven useful in many domains, including statistical physics, where it describes the behavior of gases or magnets, in computer vision for image segmentation, and in social network analysis.

2 Graphical model selection

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 (e.g., gene microarrays, social networks etc.), the model dimension $p$ is comparable or larger than the sample size $n$, so that the relevance of such “fixed $p$” asymptotics is doubtful. Accordingly, the goal of this paper is to develop the broader notion of high-dimensional consistency, in which both the model dimension and the sample size are allowed to increase, and we study the scaling conditions 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

Note that 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 $\sign(\theta^{*}_{rt})$ for all $t\in\mathcal{N}(r)$. To capture both the neighborhood structure and sign pattern, we define the signed neighborhood set

The next step is to 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^{*}_{\backslash r}$, we consider the structure of the conditional distribution of $X_{r}$ given the other variables $X_{\backslash r}=\{X_{t}\,\mid\,t\in V\backslash\{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_{\backslash r}$ play the role of the covariates.

, where $\lambda_{n}>0$ is a regularization parameter, to be specified by the user, and

is the rescaled negative log likelihood. (The rescaling factor $1/n$ in this definition is for later theoretical convenience.) Following some algebraic manipulation, the regularized negative log likelihood can be written as

Accordingly, let $\widehat{\theta}^{n}_{\backslash r}$ be an element of the minimizing set of problem (6). Although $\widehat{\theta}^{n}_{\backslash r}$ need not be unique in general since the problem (6) need not be strictly convex, our analysis shows that in the regime of interest, this minimizer $\widehat{\theta}^{n}_{\backslash r}$ is indeed unique. We use $\widehat{\theta}^{n}_{\backslash 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{G}=G(p,d)\}$, if $\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 (8), one for each node $r$ of the graph, agree with the true neighborhoods, so that the full graph $G(p,d)$ is estimated consistently. In this section, we begin by stating the assumptions that underlie our main result, 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 detail 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 exponentially fast rates, we can then use a union bound over all $p$ nodes of the graph to conclude that consistent graph selection $\{\widehat{G}=G(p,d)\}$ is also achieved.

For future reference, we calculate 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: 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., non-neighbors 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_{\operatorname{{\scriptstyle{}}}}\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 a sequence of graphs $\{G(p,d)\}$ such that conditions $A1$ and $A2$ are satisfied by the population Fisher information matrices $Q^{*}$. If the sample size $n$ satisfies

for some constant $L$, and the minimum value $\theta^{*}_{\operatorname{min}}$ decays no faster than ${\mathcal{O}}(1/d)$, then for the regularization sequence $\lambda_{n}=2\sqrt{\frac{\log p}{n}}$, the estimated graph $\widehat{G}(\lambda_{n})$ obtained by neighborhood logistic regression satisfies

For model selection in graphical models, one is typically interested in node degrees $d$ that remain bounded (e.g., $d={\mathcal{O}}(1)$), or that grow only weakly with graph size (say $d=o(p$). In such cases, the growth condition (15) allows the number of observations to be substantially smaller than the graph size, i.e., the “large $p$, small $n$” regime. In particular, the graph size $p$ can grow exponentially with the number of observations (i.e, $p(n)=\exp(n^{\alpha})$ for some $\alpha\in(0,1)$.

In terms of the choice of regularization, the sequence $\lambda_{n}$ needs to satisfy the following conditions:

Under the growth condition (15), the choice $\lambda_{n}=2\sqrt{\frac{\log p}{n}}$ suffices as long as $\theta^{*}_{\operatorname{min}}$ decays no faster than ${\mathcal{O}}(1/d)$.

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 (A1) mutual incoherence (A2) conditions hold for the sample Fisher information matrix

then the growth condition (15) and choice of $\lambda_{n}$ from Theorem 1 are sufficient to ensure that the graph is recovered with high probability. Interestingly, our analysis shows that if the conditions are imposed directly on the sample Fisher information matrices and $\theta^{*}_{\operatorname{min}}=\Theta(1)$, then the weaker growth condition $n=\Omega(d^{2}\log(p))$ suffices for asymptotically exact graph recovery.

The second part of the analysis, provided in Section 5, is devoted to showing that under the specified growth conditions (A3), 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}$. While 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$, the delicacy is that we require controlling this convergence over subsets of increasing size. The analysis therefore requires some large-deviations bounds, so as to provide exponential control on the rates of convergence.

3 Primal-dual witness for graph recovery

At a high-level, at the core of our proof lies the notion of a primal-dual witness. In particular, we explicitly construct an optimal primal-dual pair, namely, a primal solution $\widehat{\theta}$, along with an associated subgradient vector $\widehat{z}$ (which can be interpreted as a dual solution), such that the Karush-Kuhn-Tucker (KKT) conditions associated with the convex program (6) are satisfied. Moreover, we show that under the stated assumptions on $(n,p,d)$, the primal-dual pair $(\widehat{\theta},\widehat{z})$ can be constructed such that they act as a witness—that is, a certificate guaranteeing that the method correctly recovers the graph structure.

Let us write the convex program (6) in the form

is the negative log likelihood associated with the logistic regression model. The KKT conditions associated with this model can be expressed as follows

where $\alpha^{*}_{v}\geq 0$ is the Lagrange multiplier associated with the constraint ${\vec{v}}^{T}\theta\leq C$.

The KKT conditions (20) and (21) must be satisfied by any optimal pair $(\widehat{\theta},\widehat{z})$ to the convex program (18). In order for this primal-dual pair to correctly specify the graph structure, we require furthermore that the following properties are satisfied:

We now construct our witness pair $(\widehat{\theta},\widehat{z})$ as follows. First, we set $\widehat{\theta}_{S}$ as the minimizer of the partial penalized likelihood,

and set $\widehat{z}_{S}=\sign(\widehat{\theta}_{S})$. We then set $\widehat{\theta}_{{S^{c}}}=0$ so that condition (23b) holds. Finally, we obtain $\widehat{z}_{{{S^{c}}}}$ from equation (20) by plugging in the values of $\widehat{\theta}$ and $\widehat{z}_{S}$. Thus, our construction satisfies conditions (23b) and (20). The remainder of the analysis consists of showing that our conditions on $(n,p,d)$ imply that, with high-probability, the remaining conditions (23a) and (21) are satisfied.

This strategy is justified by the following lemma, which provides sufficient conditions for shared sparsity and uniqueness of the optimal solution:

By Lagrangian duality, the penalized problem (18) 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. Consider the representation of $\widehat{z}$ as the convex combination (22) of sign vectors ${\vec{v}}\in\{-1,+1\}^{p-1}$, where the weights $\alpha^{*}_{v}$ are non-negative and sum to one. Since $\alpha^{*}$ is an optimal vector of Lagrange multipliers for the optimal primal solution $\widehat{\theta}$, it follows (Bertsekas 1995) that any other optimal primal solution $\widetilde{\theta}$ must minimize the associated Lagrangian (i.e., satisfy equation (20)), and moreover must satisfy the complementary slackness conditions $\alpha^{*}_{v}[{\vec{v}}^{T}\theta-C]=0$ for all sign vectors $v$. But these conditions imply that $\widehat{z}^{T}\widetilde{\theta}=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. ∎

In our primal-dual witness proof, we exploit this lemma by constructing a primal-dual pair $(\widehat{\theta},\widehat{z})$ such that $\|\widehat{z}_{S^{c}}\|_{\infty}<1$. 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 the primal solution $\widehat{\theta}$ is guaranteed to be unique.

Analysis under sample Fisher matrix assumptions

We begin by establishing model selection consistency when assumptions are imposed directly on the sample Fisher matrix $Q^{n}$, as opposed to on the population matrix $Q^{*}$, as in Theorem 1. In particular, we define the “good event”

If $n>Ld^{2}\log(p)$ for a suitably large constant $L$, and the minimum value $\theta^{*}_{\operatorname{min}}$ decays no faster than ${\mathcal{O}}(1/\sqrt{d})$, then for the regularization sequence $\lambda_{n}=2\sqrt{\frac{\log p}{n}}$, the estimated graph $\widehat{G}(\lambda_{n})$ obtained by neighborhood logistic regression satisfies

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 A. The central object is the following expansion, obtained by re-writing the zero-subgradient condition as

where we have introduced the short-hand notation $W^{n}$ for the $(p-1)$-vector

with $\bar{\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:

If $n\lambda_{n}^{2}>\log(p)$, then for the specified mutual incoherence parameter $\alpha_{\operatorname{{\scriptstyle{}}}}\in(0,1]$, we have

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

If $\lambda_{n}d\leq\frac{C_{min}^{2}}{10D_{\operatorname{max}}}$, then as $n\rightarrow+\infty$, we have

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

Our final technical lemma provides control on the the remainder term (29):

If $n\lambda^{2}_{n}>\log(p)$ and $d\lambda_{n}$ is sufficiently small, then for mutual incoherence parameter $\alpha_{\operatorname{{\scriptstyle{}}}}\in(0,1]$, we have

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

2 Proof of Proposition 1

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

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

Correct sign recovery:

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

where we recall the notation $\theta^{*}_{\operatorname{min}}:=\min_{(r,t)\in E}|\theta^{*}_{rt}|$. From Lemma 3, we have $\|\theta_{S}-\theta^{*}_{S}\|_{2}={\mathcal{O}}_{p}(\sqrt{d}\lambda_{n})$, so that

Since $\theta^{*}_{\operatorname{min}}$ decays no faster than $\Theta(1/\sqrt{d})$, the right-hand side is upper bounded by ${\mathcal{O}}(\lambda_{n}d)$, which can be made smaller than $1$ by choosing $\lambda_{n}$ sufficiently small, as asserted in Proposition 1.

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., Davidson and Szarek 2001), since the elements of $Q^{n}$ are highly dependent.

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 (13) with parameter $\alpha_{\operatorname{{\scriptstyle{}}}}\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 begin by taking note of a simple bound to be used repeatedly throughout our arguments. By definition of the matrices $Q^{n}(\theta)$ and $Q(\theta)$ (see equations (17) and (10)), 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 (Hoeffding 1963), for any indices $j,k=1,\ldots,d$ and for any $\epsilon>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 (Horn and Johnson 1985), we have

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

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

which obeys the same upper bound (44), 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_{\operatorname{{\scriptstyle{}}}}}{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 (42), 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 B for the proof of these claims.

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

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

where we have used the incoherence assumption $A2$. Using the bound (41b) 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\left(-Kn/d^{2}+2\log(d)\right)$. Next, applying the bound (46b) with $\delta=c/\sqrt{d}$, we conclude that with probability greater than $1-2\exp\left(-Knc^{2}/d^{3}+\log(d)\right)$, we have

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

Control of second term:

We then apply bound (46a) with $\delta=\frac{\alpha_{\operatorname{{\scriptstyle{}}}}}{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 (46a) and (46b), both with $\delta=\sqrt{\alpha_{\operatorname{{\scriptstyle{}}}}/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

We performed experiments for three different classes of graphs: four-nearest neighbor lattices, (b) eight-nearest neighbor lattices, and (c) star-shaped graphs, as illustrated in Figure 1.

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, 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\}$, as well as both mixed and attractive couplings. Notice how once again the curves for different problem sizes are all well-aligned, consistent with the prediction of Theorem 1.

For our last set of experiments, we investigated the performance of our method for a class of graphs with unbounded maximum degree $d$. In particular, we constructed star-shaped graphs with $p$ vertices by designating one node as the spoke, and connecting it to $d<(p-1)$ of its neighbors. For linear sparsity, we chose $d=\lceil 0.1p\rceil$, whereas for logarithmic sparsity we choose $d=\lceil\log(p)\rceil$. We again studied 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.

Conclusion

Research supported in part by NSF grants IIS-0427206 and CCF-0625879 (PR and JL), NSF grants DMS-0605165 and CCF-0545862 (PR and MJW), and a Siebel Scholarship (PR).

Appendix A Proofs for Section 4.1

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

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

for some constant $K$. Finally, applying a union bound over the indices $u$ of $W^{n}$ yields

A.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}$ with probability converging to one from Lemma 2.

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

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 $|u_{S}^{T}x^{(i)}_{S}|\leq\sqrt{d}\|u_{S}\|_{2}=M\lambda_{n}d$. Moreover, we have $\|\frac{1}{n}\sum_{i=1}^{n}\left(x^{(i)}_{S})^{T}y\right)\leq\|\frac{1}{n}\sum_{i=1}^{n}x^{(i)}_{S}(x^{(i)}_{S})^{T}\|_{2}\;\leq\;D_{\operatorname{max}}$ by assumption. Combining these pieces, we obtain

where the last inequality follows as long as $\lambda_{n}d\leq\frac{C_{min}}{2D_{\operatorname{max}}M}$. We have thus shown that

with probability converging to one, as long as $\lambda_{n}d$ is sufficiently small.

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}$. Moreover, for this choice of $M$, we have that $\lambda_{n}d$ must be upper bounded by $\frac{C_{min}}{2D_{\operatorname{max}}M}=\frac{C_{min}^{2}}{10D_{\operatorname{max}}}$, as assumed in the lemma statement.

A.3 Proof of Lemma 4

We first show that the remainder term $R^{n}$ satisfies the bound $\|R^{n}\|_{\infty}\leq D_{\operatorname{max}}\|\widehat{\theta}_{S}-\theta^{*}_{S}\|_{2}^{2}$. Then the result of Lemma 3—namely, that $\|\widehat{\theta}_{S}-\theta^{*}_{S}\|_{2}=\mathcal{O}_{p}(\lambda_{n}\sqrt{d})$—can be used to conclude that $\frac{\|R^{n}\|_{\infty}}{\lambda_{n}}=\mathcal{O}_{p}(\lambda_{n}d)$, which suffices to guarantee the claim of Lemma 4.

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

for some point $\bar{\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\backslash r}\theta_{rt}x_{t})$. By the chain rule and another application of the mean value theorem, we then have

where $\bar{\bar{\theta}}^{(j)}$ is another point on the line joining $\widehat{\theta}$ and ${{\theta^{*}}}$. Setting $a_{i}:=\{g^{\prime}\left(\bar{\bar{\theta}}^{(j)T}x^{(i)}\right)x^{(i)}_{j}\}$ and $b_{i}:=\{[\bar{\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 B 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, we have

from which the claim (46a) follows by setting $\epsilon=\delta/d$ in the Hoeffding bound (43). The proof of bound (46b) is analogous, with the pre-factor $(p-d)$ replaced by $d$.

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

for a constants $B$. Moreover, from equation (44), we have

References