Learning the Structure of Mixed Graphical Models

# Learning Mixed Graphical Models

## Abstract

We consider the problem of learning the structure of a pairwise graphical model over continuous and discrete variables. We present a new pairwise model for graphical models with both continuous and discrete variables that is amenable to structure learning. In previous work, authors have considered structure learning of Gaussian graphical models and structure learning of discrete models. Our approach is a natural generalization of these two lines of work to the mixed case. The penalization scheme involves a novel symmetric use of the group-lasso norm and follows naturally from a particular parametrization of the model.

## 1 Introduction

Many authors have considered the problem of learning the edge structure and parameters of sparse undirected graphical models. We will focus on using the regularizer to promote sparsity. This line of work has taken two separate paths: one for learning continuous valued data and one for learning discrete valued data. However, typical data sources contain both continuous and discrete variables: population survey data, genomics data, url-click pairs etc. For genomics data, in addition to the gene expression values, we have attributes attached to each sample such as gender, age, ethniticy etc. In this work, we consider learning mixed models with both continuous variables and discrete variables.

For only continuous variables, previous work assumes a multivariate Gaussian (Gaussian graphical) model with mean and inverse covariance . is then estimated via the graphical lasso by minimizing the regularized negative log-likelihood . Several efficient methods for solving this can be found in Friedman et al. (2008a); Banerjee et al. (2008). Because the graphical lasso problem is computationally challenging, several authors considered methods related to the pseudolikelihood (PL) and nodewise regression (Meinshausen and Bühlmann, 2006; Friedman et al., 2010; Peng et al., 2009). For discrete models, previous work focuses on estimating a pairwise Markov random field of the form . The maximum likelihood problem is intractable for models with a moderate to large number of variables (high-dimensional) because it requires evaluating the partition function and its derivatives. Again previous work has focused on the pseudolikelihood approach (Guo et al., 2010; Schmidt, 2010; Schmidt et al., 2008; Höfling and Tibshirani, 2009; Jalali et al., 2011; Lee et al., 2006; Ravikumar et al., 2010).

Our main contribution here is to propose a model that connects the discrete and continuous models previously discussed. The conditional distributions of this model are two widely adopted and well understood models: multiclass logistic regression and Gaussian linear regression. In addition, in the case of only discrete variables, our model is a pairwise Markov random field; in the case of only continuous variables, it is a Gaussian graphical model. Our proposed model leads to a natural scheme for structure learning that generalizes the graphical Lasso. Here the parameters occur as singletons, vectors or blocks, which we penalize using group-lasso norms, in a way that respects the symmetry in the model. Since each parameter block is of different size, we also derive a calibrated weighting scheme to penalize each edge fairly. We also discuss a conditional model (conditional random field) that allows the output variables to be mixed, which can be viewed as a multivariate response regression with mixed output variables. Similar ideas have been used to learn the covariance structure in multivariate response regression with continuous output variables Witten and Tibshirani (2009); Kim et al. (2009); Rothman et al. (2010).

In Section 2, we introduce our new mixed graphical model and discuss previous approaches to modeling mixed data. Section 3 discusses the pseudolikelihood approach to parameter estimation and connections to generalized linear models. Section 4 discusses a natural method to perform structure learning in the mixed model. Section 5 presents the calibrated regularization scheme, Section 6 discusses the consistency of the estimation procedures, and Section 7 discusses two methods for solving the optimization problem. Finally, Section 8 discusses a conditional random field extension and Section 9 presents empirical results on a census population survey dataset and synthetic experiments.

## 2 Mixed Graphical Model

We propose a pairwise graphical model on continuous and discrete variables. The model is a pairwise Markov random field with density proportional to

 exp(p∑s=1p∑t=1−12βstxsxt+p∑s=1αsxs+p∑s=1q∑j=1ρsj(yj)xs+q∑j=1q∑r=1ϕrj(yr,yj)). (1)

Here denotes the th of continuous variables, and the th of discrete variables. The joint model is parametrized by 3. The discrete takes on states. The model parameters are continuous-continuous edge potential, continuous node potential, continuous-discrete edge potential, and discrete-discrete edge potential.

The two most important features of this model are:

1. the conditional distributions are given by Gaussian linear regression and multiclass logistic regressions;

2. the model simplifies to a multivariate Gaussian in the case of only continuous variables and simplifies to the usual discrete pairwise Markov random field in the case of only discrete variables.

The conditional distributions of a graphical model are of critical importance. The absence of an edge corresponds to two variables being conditionally independent. The conditional independence can be read off from the conditional distribution of a variable on all others. For example in the multivariate Gaussian model, is conditionally independent of iff the partial correlation coefficient is . The partial correlation coefficient is also the regression coefficient of in the linear regression of on all other variables. Thus the conditional independence structure is captured by the conditional distributions via the regression coefficient of a variable on all others. Our mixed model has the desirable property that the two type of conditional distributions are simple Gaussian linear regressions and multiclass logistic regressions. This follows from the pairwise property in the joint distribution. In more detail:

1. The conditional distribution of given the rest is multinomial, with probabilities defined by a multiclass logistic regression where the covariates are the other variables and (denoted collectively by in the right-hand side):

 p(yr=k|y∖r,x;Θ)=exp(ωTkz)Lr∑l=1exp(ωTlz)=exp(ω0k+∑jωkjzj)∑Lrl=1exp(ω0l+∑jωljzj) (2)

Here we use a simplified notation, which we make explicit in Section 3.1. The discrete variables are represented as dummy variables for each state, e.g. , and for continuous variables .

2. The conditional distribution of given the rest is Gaussian, with a mean function defined by a linear regression with predictors and .

 E(xs|x∖s,yr;Θ) =ωTz=ω0+∑jzjωj (3) p(xs|x∖s,yr;Θ) =1√2πσsexp(−12σ2s(xs−ωTz)2).

As before, the discrete variables are represented as dummy variables for each state and for continuous variables .

The exact form of the conditional distributions (2) and (3) are given in (11) and (10) in Section 3.1, where the regression parameters are defined in terms of the parameters .

The second important aspect of the mixed model is the two special cases of only continuous and only discrete variables.

1. Continuous variables only. The pairwise mixed model reduces to the familiar multivariate Gaussian parametrized by the symmetric positive-definite inverse covariance matrix and mean ,

 p(x)∝exp(−12(x−B−1α)TB(x−B−1α)).
2. Discrete variables only. The pairwise mixed model reduces to a pairwise discrete (second-order interaction) Markov random field,

 Missing or unrecognized delimiter for \right

Although these are the most important aspects, we can characterize the joint distribution further. The conditional distribution of the continuous variables given the discrete follow a multivariate Gaussian distribution, . Each of these Gaussian distributions share the same inverse covariance matrix but differ in the mean parameter, since all the parameters are pairwise. By standard multivariate Gaussian calculations,

 p(x|y) =N(B−1γ(y),B−1) (4) {γ(y)}s =αs+∑jρsj(yj) (5) p(y) ∝exp(q∑j=1j∑r=1ϕrj(yr,yj)+12γ(y)TB−1γ(y)) (6)

Thus we see that the continuous variables conditioned on the discrete are multivariate Gaussian with common covariance, but with means that depend on the value of the discrete variables. The means depend additively on the values of the discrete variables since . The marginal has a known form, so for models with few number of discrete variables we can sample efficiently.

### 2.1 Related work on mixed graphical models

Lauritzen (1996) proposed a type of mixed graphical model, with the property that conditioned on discrete variables, . The homogeneous mixed graphical model enforces common covariance, . Thus our proposed model is a special case of Lauritzen’s mixed model with the following assumptions: common covariance, additive mean assumptions and the marginal factorizes as a pairwise discrete Markov random field. With these three assumptions, the full model simplifies to the mixed pairwise model presented. Although the full model is more general, the number of parameters scales exponentially with the number of discrete variables, and the conditional distributions are not as convenient. For each state of the discrete variables there is a mean and covariance. Consider an example with binary variables and continuous variables; the full model requires estimates of mean vectors and covariance matrices in dimensions. Even if the homogeneous constraint is imposed on Lauritzen’s model, there are still mean vectors for the case of binary discrete variables. The full mixed model is very complex and cannot be easily estimated from data without some additional assumptions. In comparison, the mixed pairwise model has number of parameters and allows for a natural regularization scheme which makes it appropriate for high dimensional data.

An alternative to the regularization approach that we take in this paper, is the limited-order correlation hypothesis testing method Tur and Castelo (2012). The authors develop a hypothesis test via likelihood ratios for conditional independence. However, they restrict to the case where the discrete variables are marginally independent so the maximum likelihood estimates are well-defined for .

There is a line of work regarding parameter estimation in undirected mixed models that are decomposable: any path between two discrete variables cannot contain only continuous variables. These models allow for fast exact maximum likelihood estimation through node-wise regressions, but are only applicable when the structure is known and (Edwards, 2000). There is also related work on parameter learning in directed mixed graphical models. Since our primary goal is to learn the graph structure, we forgo exact parameter estimation and use the pseudolikelihood. Similar to the exact maximum likelihood in decomposable models, the pseudolikelihood can be interpreted as node-wise regressions that enforce symmetry.

To our knowledge, this work is the first to consider convex optimization procedures for learning the edge structure in mixed graphical models.

## 3 Parameter Estimation: Maximum Likelihood and Pseudolikelihood

Given samples , we want to find the maximum likelihood estimate of . This can be done by minimizing the negative log-likelihood of the samples:

 ℓ(Θ) =−n∑i=1logp(xi,yi;Θ) where (7) logp(x,y;Θ) =p∑s=1p∑t=1−12βstxsxt+p∑s=1αsxs+p∑s=1q∑j=1ρsj(yj)xs +q∑j=1j∑r=1ϕrj(yr,yj)−logZ(Θ) (8)

The negative log-likelihood is convex, so standard gradient-descent algorithms can be used for computing the maximum likelihood estimates. The major obstacle here is , which involves a high-dimensional integral. Since the pairwise mixed model includes both the discrete and continuous models as special cases, maximum likelihood estimation is at least as difficult as the two special cases, the first of which is a well-known computationally intractable problem. We defer the discussion of maximum likelihood estimation to Appendix .5.

### 3.1 Pseudolikelihood

The pseudolikelihood method Besag (1975) is a computationally efficient and consistent estimator formed by products of all the conditional distributions:

 ~ℓ(Θ|x,y)=−p∑s=1logp(xs|x∖s,y;Θ)−q∑r=1logp(yr|x,y∖r;Θ) (9)

The conditional distributions and take on the familiar form of linear Gaussian and (multiclass) logistic regression, as we pointed out in (2) and (3). Here are the details:

• The conditional distribution of a continuous variable is Gaussian with a linear regression model for the mean, and unknown variance.

 p(xs|x∖s,y;Θ)=√βss√2πexp⎛⎝−βss2(αs+∑jρsj(yj)−∑t≠sβstxtβss−xs)2⎞⎠ (10)
• The conditional distribution of a discrete variable with states is a multinomial distribution, as used in (multiclass) logistic regression. Whenever a discrete variable is a predictor, each of its levels contribute an additive effect; continuous variables contribute linear effects.

 p(yr|y∖r,,x;Θ)=exp(∑sρsr(yr)xs+ϕrr(yr,yr)+∑j≠rϕrj(yr,yj))∑Lrl=1exp(∑sρsr(l)xs+ϕrr(l,l)+∑j≠rϕrj(l,yj)) (11)

Taking the negative log of both gives us

 −logp(xs|x∖s,y;Θ) =−12logβss+βss2⎛⎝αsβss+∑jρsj(yj)βss−∑t≠sβstβssxt−xs⎞⎠2 (12) −logp(yr|y∖r,,x;Θ) =−logexp(∑sρsr(yr)xs+ϕrr(yr,yr)+∑j≠rϕrj(yr,yj))∑Lrl=1exp(∑sρsr(l)xs+ϕrr(l,l)+∑j≠rϕrj(l,yj)) (13)

A generic parameter block, , corresponding to an edge appears twice in the pseudolikelihood, once for each of the conditional distributions and .

###### Proposition 1.

The negative log pseudolikelihood in (9) is jointly convex in all the parameters over the region .

We prove Proposition 1 in Appendix .3.

### 3.2 Separate node-wise regression

A simple approach to parameter estimation is via separate node-wise regressions; a generalized linear model is used to estimate for each . Separate regressions were used in Meinshausen and Bühlmann (2006) for the Gaussian graphical model and Ravikumar et al. (2010) for the Ising model. The method can be thought of as an asymmetric form of the pseudolikelihood since the pseudolikelihood enforces that the parameters are shared across the conditionals. Thus the number of parameters estimated in the separate regression is approximately double that of the pseudolikelihood, so we expect that the pseudolikelihood outperforms at low sample sizes and low regularization regimes. The node-wise regression was used as our baseline method since it is straightforward to extend it to the mixed model. As we predicted, the pseudolikelihood or joint procedure outperforms separate regressions; see top left box of Figures 4 and 5. Liu and Ihler (2012, 2011) confirm that the separate regressions are outperformed by pseudolikelihood in numerous synthetic settings.

Concurrent work of Yang et al. (2012, 2013) extend the separate node-wise regression model from the special cases of Gaussian and categorical regressions to generalized linear models, where the univariate conditional distribution of each node is specified by a generalized linear model (e.g. Poisson, categorical, Gaussian). By specifying the conditional distributions, Besag (1974) show that the joint distribution is also specified. Thus another way to justify our mixed model is to define the conditionals of a continuous variable as Gaussian linear regression and the conditionals of a categorical variable as multiple logistic regression and use the results in Besag (1974) to arrive at the joint distribution in (1). However, the neighborhood selection algorithm in Yang et al. (2012, 2013) is restricted to models of the form In particular, this procedure cannot be applied to edge selection in our pairwise mixed model in (1) or the categorical model in (2) with greater than 2 states. Our baseline method of separate regressions is closely related to the neighborhood selection algorithm they proposed; the baseline can be considered as a generalization of Yang et al. (2012, 2013) to allow for more general pairwise interactions with the appropriate regularization to select edges. Unfortunately, the theoretical results in Yang et al. (2012, 2013) do not apply to the baseline nodewise regression method, nor the joint pseudolikelihood.

## 4 Conditional Independence and Penalty Terms

In this section, we show how to incorporate edge selection into the maximum likelihood or pseudolikelihood procedures. In the graphical representation of probability distributions, the absence of an edge corresponds to a conditional independency statement that variables and are conditionally independent given all other variables (Koller and Friedman, 2009). We would like to maximize the likelihood subject to a penalization on the number of edges since this results in a sparse graphical model. In the pairwise mixed model, there are 3 type of edges

1. is a scalar that corresponds to an edge from to . implies and are conditionally independent given all other variables. This parameter is in two conditional distributions, corresponding to either or is the response variable, and .

2. is a vector of length . If for all values of , then and are conditionally independent given all other variables. This parameter is in two conditional distributions, corresponding to either or being the response variable: and .

3. is a matrix of size . If for all values of and , then and are conditionally independent given all other variables. This parameter is in two conditional distributions, corresponding to either or being the response variable, and .

For edges that involve discrete variables, the absence of that edge requires that the entire matrix or vector is . The form of the pairwise mixed model motivates the following regularized optimization problem

 minimizeΘ ℓλ(Θ)=ℓ(Θ)+λ(∑s

All parameters that correspond to the same edge are grouped in the same indicator function. This problem is non-convex, so we replace the sparsity and group sparsity penalties with the appropriate convex relaxations. For scalars, we use the absolute value ( norm), for vectors we use the norm, and for matrices we use the Frobenius norm. This choice corresponds to the standard relaxation from group to group (group lasso) norm (Bach et al., 2011; Yuan and Lin, 2006).

 Missing or unrecognized delimiter for \left (15)

## 5 Calibrated regularizers

In (15) each of the group penalties are treated as equals, irrespective of the size of the group. We suggest a calibration or weighting scheme to balance the load in a more equitable way. We introduce weights for each group of parameters and show how to choose the weights such that each parameter set is treated equally under , the fully-factorized independence model 4

 minimizeΘ ℓ(Θ)+λ(p∑t=1s−1∑t=1wst|βst|+p∑s=1q∑j=1wsj∥∥ρsj∥∥2+q∑j=1j−1∑r=1wrj∥∥ϕrj∥∥F) (16)

Based on the KKT conditions (Friedman et al., 2007), the parameter group is non-zero if

 ∥∥∥∂ℓ∂θg∥∥∥>λwg

where and represents one of the parameter groups and its corresponding weight. Now can be viewed as a generalized residual, and for different groups these are different dimensions—e.g. scalar/vector/matrix. So even under the independence model (when all terms should be zero), one might expect some terms to have a better random chance of being non-zero (for example, those of bigger dimensions). Thus for all parameters to be on equal footing, we would like to choose the weights such that

 EpF∥∥∥∂ℓ∂θg∥∥∥=constant% ×wg

However, it is simpler to compute in closed form , so we choose

 wg∝ ⎷EpF∥∥∥∂ℓ∂θg∥∥∥2

where is the fully factorized (independence) model. In Appendix .6, we show that the weights can be chosen as

 wst =σsσt wsj =σs√∑apa(1−pa) wrj =√∑apa(1−pa)∑bqb(1−qb)

is the standard deviation of the continuous variable . and . For all types of parameters, the weight has the form of , where represents a generic variable and is the variance-covariance matrix of .

## 6 Model Selection Consistency

In this section, we study the model selection consistency, the correct edge set is selected and the parameter estimates are close to the truth, of pseudolikelihood and maximum likelihood. We will see that the consistency can be established using the framework first developed in Ravikumar et al. (2010) and later extended to general m-estimators by Lee et al. (2013). The proofs in this section are omitted since they follow from a straightforward application of the results in Lee et al. (2013); the results are stated for the mixed model to show that under certain conditions the estimation procedures are model selection consistent. We also only consider the uncalibrated regularizers to simplify the notation, but it is straightforward to adapt to the calibrated regularizer case.

First, we define some notation. Recall that is the vector of parameters being estimated , be the true parameters that estimated the model, and . Both estimation procedures can be written as a convex optimization problem of the form

 minimize ℓ(Θ)+λ∑g∈G∥∥Θg∥∥2 (17)

where is one of the two log-likelihoods. The regularizer

 ∑g∈G∥∥Θg∥∥=λ(p∑s=1s−1∑t=1|βst|+p∑s=1q∑j=1∥∥ρsj∥∥2+q∑j=1j−1∑r=1∥∥ϕrj∥∥F).

The set indexes the edges , , and , and is one of the three types of edges.

It is difficult to establish consistency results for the problem in Equation (17) because the parameters are non-identifiable. This is because is constant with respect to the change of variables and similarly for , so we cannot hope to recover . A popular fix for this issue is to drop the last level of and , so they are only indicators over levels instead of levels. This allows for the model to be identifiable, but it results in an asymmetric formulation that treats the last level differently from other levels. Instead, we will maintain the symmetric formulation by introducing constraints. Consider the problem

 minimizeΘ ℓ(Θ)+λ∑g∈G∥∥Θg∥∥2 (18) subject to CΘ=0.

The matrix constrains the optimization variables such that

 ∑yjρsj(yj)=0 ∑yjϕrj(yr,yj)=0.

The group regularizer implicitly enforces the same set of constraints, so the optimization problems of Equation (18) and Equation (17) have the same solutions. For our theoretical results, we will use the constrained formulation of Equation (18), since it is identifiable.

We first state some definitions and two assumptions from Lee et al. (2013) that are necessary to present the model selection consistency results. Let and represent the active and inactive groups in , so for any and for any . The sets associated with the active and inactive groups are defined as

 A ={Θ∈Rd:maxg∈G∥∥Θg∥∥2≤1and ∥∥Θg∥∥2=0,g∈I} I ={Θ∈Rd:maxg∈G∥∥Θg∥∥2≤1and ∥∥Θg∥∥2=0,g∈A}.

Let and be the orthogonal projector onto the subspace . The two assumptions are

1. Restricted Strong Convexity. We assume that

 supv∈MvT∇2ℓ(Θ)vvTv≥m (19)

for all . Since is lipschitz continuous, the existence of a constant that satisfies (19) is implied by the pointwise restricted convexity

 supv∈MvT∇2ℓ(Θ⋆)vvTv≥~m.

For convenience, we will use the former.

2. Irrepresentable condition. There exist such that

 supz∈AV(PM⊥(∇2ℓ(Θ⋆)PM(PM∇2ℓ(Θ⋆)PM)†PMz−z))<1−τ, (20)

where is the infimal convolution of , the gauge of set , and :

 V(z)=infz=u1+u2{ρI(u1)+\mathbbm1Null(C)⊥(u2)}.

Restricted strong convexity is a standard assumption that ensures the parameter is uniquely determined by the value of the likelihood function. Without this, there is no hope of accurately estimating . It is only stated over a subspace which can be much smaller than . The Irrepresentable condition is a more stringent condition. Intuitively, it requires that the active variables not be overly dependent on the inactive variables. Although the exact form of the condition is not enlightening, it is known to be ”almost” necessary for model selection consistency in the lasso (Zhao and Yu, 2006) and a common assumption in other works that establish model selection consistency (Ravikumar et al., 2010; Jalali et al., 2011; Peng et al., 2009). We also define the constants that appear in the theorem:

1. Lipschitz constants and . Let be the log-partition function. and are twice continuously differentiable functions, so their gradient and hessian are locally Lipschitz continuous in a ball of radius around :

 Missing or unrecognized delimiter for \right ∥∥∇2ℓ(Θ1)−∇2ℓ(Θ2)∥∥2≤L2∥Θ1−Θ2∥2, Θ1,Θ2∈Br(Θ⋆)
2. Let satisfy

 supz∈A∪IV(PM⊥(∇2ℓ(Θ⋆)PM(PM∇2ℓ(Θ⋆)PM)†PMz−z))<¯τ.

is a continuous function of , so a finite exists.

###### Theorem 2.

Suppose we are given samples from the mixed model with unknown parameters . If we select

 λ=2√256L1¯ττ√(maxg∈G|g|)log|G|n

and the sample size is larger than

 Missing or unrecognized delimiter for \left

then, with probability at least , the optimal solution to (17) is unique and model selection consistent,

1. and .

###### Remark 3.

The same theorem applies to both the maximum likelihood and pseudolikelihood estimators. For the maximum likelihood, the constants can be tightened; everywhere appears can be replaced by and the theorem remains true. However, the values of are different for the two methods. For the maximum likelihood, the gradient of the log-partition and hessian of the log-likelihood do not depend on the samples. Thus the constants are completely determined by and the likelihood. For the pseudolikelihood, the values of depend on the samples, and the theorem only applies if the assumptions are made on sample quantities; thus, the theorem is less useful in practice when applied to the pseudolikelihood. This is similar to the situation in Yang et al. (2013), where assumptions are made on sample quantities.

## 7 Optimization Algorithms

In this section, we discuss two algorithms for solving (15): the proximal gradient and the proximal newton methods. This is a convex optimization problem that decomposes into the form , where is smooth and convex and is convex but possibly non-smooth. In our case is the negative log-likelihood or negative log-pseudolikelihood and are the group sparsity penalties.

Block coordinate descent is a frequently used method when the non-smooth function is the or group . It is especially easy to apply when the function is quadratic, since each block coordinate update can be solved in closed form for many different non-smooth (Friedman et al., 2007). The smooth in our particular case is not quadratic, so each block update cannot be solved in closed form. However in certain problems (sparse inverse covariance), the update can be approximately solved by using an appropriate inner optimization routine (Friedman et al., 2008b).

Problems of this form are well-suited for the proximal gradient and accelerated proximal gradient algorithms as long as the proximal operator of can be computed (Combettes and Pesquet, 2011; Beck and Teboulle, 2010)

 Missing or unrecognized delimiter for \right (21)

For the sum of group sparsity penalties considered, the proximal operator takes the familiar form of soft-thresholding and group soft-thresholding (Bach et al., 2011). Since the groups are non-overlapping, the proximal operator simplifies to scalar soft-thresholding for and group soft-thresholding for and .

The class of proximal gradient and accelerated proximal gradient algorithms is directly applicable to our problem. These algorithms work by solving a first-order model at the current iterate

 argminu f(xk)+∇f(xk)T(u−xk)+12t∥u−xk∥2+g(u) (22) =argminu 12t∥u−(xk−t∇f(xk))∥2+g(u) (23) =proxt(xk−t∇f(xk)) (24)

The proximal gradient iteration is given by where is determined by line search. The theoretical convergence rates and properties of the proximal gradient algorithm and its accelerated variants are well-established (Beck and Teboulle, 2010). The accelerated proximal gradient method achieves linear convergence rate of when the objective is strongly convex and the sublinear rate for non-strongly convex problems.

The TFOCS framework (Becker et al., 2011) is a package that allows us to experiment with different variants of the accelerated proximal gradient algorithm. The TFOCS authors found that the Auslender-Teboulle algorithm exhibited less oscillatory behavior, and proximal gradient experiments in the next section were done using the Auslender-Teboulle implementation in TFOCS.

### 7.2 Proximal Newton Algorithms

This section borrows heavily from Schmidt (2010), Schmidt et al. (2011) and Lee et al. (2012). The class of proximal Newton algorithms is a 2nd order analog of the proximal gradient algorithms with a quadratic convergence rate (Lee et al., 2012). It attempts to incorporate 2nd order information about the smooth function into the model function. At each iteration, it minimizes a quadratic model centered at

 argminu f(xk)+∇f(xk)T(u−xk)+12t(u−xk)TH(u−xk)+g(u) (25) =argminu 12t(u−xk+tH−1∇f(xk))TH(u−xk+tH−1∇f(xk))+g(u) (26) =argminu 12t∥∥u−(xk−tH−1∇f(xk))∥∥2H+g(u) (27) :=Hproxt(xk−tH−1∇f(xk)) where % H=∇2f(xk) (28)

The operator is analogous to the proximal operator, but in the -norm. It simplifies to the proximal operator if , but in the general case of positive definite there is no closed-form solution for many common non-smooth (including and group ). However if the proximal operator of is available, each of these sub-problems can be solved efficiently with proximal gradient. In the case of separable , coordinate descent is also applicable. Fast methods for solving the subproblem include coordinate descent methods, proximal gradient methods, or Barzilai-Borwein (Friedman et al., 2007; Combettes and Pesquet, 2011; Beck and Teboulle, 2010; Wright et al., 2009). The proximal Newton framework allows us to bootstrap many previously developed solvers to the case of arbitrary loss function .

Theoretical analysis in Lee et al. (2012) suggests that proximal Newton methods generally require fewer outer iterations (evaluations of ) than first-order methods while providing higher accuracy because they incorporate 2nd order information. We have confirmed empirically that the proximal Newton methods are faster when is very large or the gradient is expensive to compute (e.g. maximum likelihood estimation). Since the objective is quadratic, coordinate descent is also applicable to the subproblems. The hessian matrix can be replaced by a quasi-newton approximation such as BFGS/L-BFGS/SR1. In our implementation, we use the PNOPT implementation (Lee et al., 2012).

### 7.3 Path Algorithm

Frequently in machine learning and statistics, the regularization parameter is heavily dependent on the dataset. is generally chosen via cross-validation or holdout set performance, so it is convenient to provide solutions over an interval of . We start the algorithm at and solve, using the previous solution as warm start, for . We find that this reduces the cost of fitting an entire path of solutions (See Figure 3). can be chosen as the smallest value such that all parameters are by using the KKT equations (Friedman et al., 2007).

## 8 Conditional Model

We can generalize our mixed model, when there are additional features , to a class of conditional random fields. Conditional models only model the conditional distribution , as opposed to the joint distribution , where are the variables of interest to the prediction task and are features.

In addition to observing and , we observe features and we build a graphical model for the conditional distribution . Consider a full pairwise model of the form (1). We then choose to only model the joint distribution over only the variables and to give us which is of the form

 p(x,y|f;Θ)= 1Z(Θ|f)exp(p∑s=1p∑t=1−12βstxsxt+p∑s=1αsxs+p∑s=1q∑j=1ρsj(yj)xs +q∑j=1j∑r=1ϕrj(yr,yj)+F∑l=1p∑s=1γlsxsfl+F∑l=1q∑r=1ηlr(yr)fl) (29)

We can also consider a more general model where each pairwise edge potential depends on the features

 p(x,y|f;Θ)=1Z(Θ|f) exp(p∑s=1p∑t=1−12βst(f)xsxt+p∑s=1αs(f)xs +p∑s=1q∑j=1ρsj(yj,f)xs+q∑j=1j∑r=1ϕrj(yr,yj,f)p∑s=1p∑t=1) (30)

(29) is a special case of this where only the node potentials depend on features and the pairwise potentials are independent of feature values. The specific parametrized form we consider is for , , and . The node potentials depend linearly on the feature values, , and .

## 9 Experimental Results

We present experimental results on synthetic data, survey data and on a conditional model.

### 9.1 Synthetic Experiments

In the synthetic experiment, the training points are sampled from a true model with continuous variables and binary variables. The edge structure is shown in Figure a. is chosen as as suggested by the theoretical results in Section 6. We see from the experimental results that recovery of the correct edge set undergoes a sharp phase transition, as expected. With samples, the pseudolikelihood is recovering the correct edge set with probability nearly . The phase transition experiments were done using the proximal Newton algorithm discussed in Section 7.2.

### 9.2 Survey Experiments

The census survey dataset we consider consists of variables, of which are continuous and are discrete: age (continuous), log-wage (continuous), year( states), sex( states),marital status ( states), race( states), education level ( states), geographic region( states), job class ( states), health ( states), and health insurance ( states). The dataset was assembled by Steve Miller of OpenBI.com from the March 2011 Supplement to Current Population Survey data. All the evaluations are done using a holdout test set of size for the survey experiments. The regularization parameter is varied over the interval at points equispaced on log-scale for all experiments.

#### Model Selection

In Figure 3, we study the model selection performance of learning a graphical model over the variables under different training samples sizes. We see that as the sample size increases, the optimal model is increasingly dense, and less regularization is needed.

#### Comparing against Separate Regressions

A sensible baseline method to compare against is a separate regression algorithm. This algorithm fits a linear Gaussian or (multiclass) logistic regression of each variable conditioned on the rest. We can evaluate the performance of the pseudolikelihood by evaluating for linear regression and for (multiclass) logistic regression. Since regression is directly optimizing this loss function, it is expected to do better. The pseudolikelihood objective is similar, but has half the number of parameters as the separate regressions since the coefficients are shared between two of the conditional likelihoods. From Figures 4 and 5, we can see that the pseudolikelihood performs very similarly to the separate regressions and sometimes even outperforms regression. The benefit of the pseudolikelihood is that we have learned parameters of the joint distribution and not just of the conditionals . On the test dataset, we can compute quantities such as conditionals over arbitrary sets of variables and marginals (Koller and Friedman, 2009). This would not be possible using the separate regressions.

#### Conditional Model

Using the conditional model (29), we model only the variables logwage, education() and jobclass(). The other variables are only used as features. The conditional model is then trained using the pseudolikelihood. We compare against the generative model that learns a joint distribution on all variables. From Figure 6, we see that the conditional model outperforms the generative model, except at small sample sizes. This is expected since the conditional distribution models less variables. At very small sample sizes and small , the generative model outperforms the conditional model. This is likely because generative models converge faster (with less samples) than discriminative models to its optimum.

#### Maximum Likelihood vs Pseudolikelihood

The maximum likelihood estimates are computable for very small models such as the conditional model previously studied. The pseudolikelihood was originally motivated as an approximation to the likelihood that is computationally tractable. We compare the maximum likelihood and maximum pseudolikelihood on two different evaluation criteria: the negative log likelihood and negative log pseudolikelihood. In Figure 7, we find that the pseudolikelihood outperforms maximum likelihood under both the negative log likelihood and negative log pseudolikelihood. We would expect that the pseudolikelihood trained model does better on the pseudolikelihood evaluation and maximum likelihood trained model does better on the likelihood evaluation. However, we found that the pseudolikelihood trained model outperformed the maximum likelihood trained model on both evaluation criteria. Although asymptotic theory suggests that maximum likelihood is more efficient than the pseudolikelihood, this analysis is applicable because of the finite sample regime and misspecified model. See Liang and Jordan (2008) for asymptotic analysis of pseudolikelihood and maximum likelihood under a well-specified model. We also observed the pseudolikelihood slightly outperforming the maximum likelihood in the synthetic experiment of Figure b.

## Acknowledgements

We would like to thank Percy Liang and Rahul Mazumder for helpful discussions. The work on consistency follows from a collaboration with Yuekai Sun and Jonathan Taylor. Jason Lee is supported by the Department of Defense (DoD) through the National Defense Science & Engineering Graduate Fellowship (NDSEG) Program, National Science Foundation Graduate Research Fellowship Program, and the Stanford Graduate Fellowship. Trevor Hastie was partially supported by grant DMS-1007719 from the National Science Foundation, and grant RO1-EB001988-15 from the National Institutes of Health.

## Appendix

### .3 Proof of Convexity

Proposition 1.The negative log pseudolikelihood in (9) is jointly convex in all the parameters over the region .

###### Proof.

To verify the convexity of , it suffices to check that each term is convex. is jointly convex in and since it is a multiclass logistic regression. We now check that is convex. is a convex function. To establish that is convex, we use the fact that is convex. Let , , and . Notice that , , , and are fixed quantities and