Sparse Gaussian Conditional Random Fields:
Algorithms, Theory, and Application to Energy Forecasting

Matt Wytock
J. Zico Kolter
Carnegie Mellon University, Pittsburgh PA

Abstract
This paper considers the sparse Gaussianconditional random field, a discriminativeextension of sparse inverse covariance estimation, where we use convex methods tolearn a high-dimensional conditional distribution of outputs given inputs. The model
has been proposed by multiple researchers
within the past year, yet previous papers
have been substantially limited in their analysis of the method and in the ability to solve
large-scale problems. In this paper, we makethree contributions: 1) we develop a secondorder active-set method which is several orders of magnitude faster than previously proposed optimization approaches for this problem, 2) we analyze the model from a theoretical standpoint, improving upon past boundswith convergence rates that depend logarithmically on the data dimension, and 3) we apply the method to large-scale energy forecasting problems, demonstrating state-of-the-art
performance on two real-world tasks.

1. Introduction
Sparse inverse covariance estimation using `1 methods(Banerjee et al., 2008), also known as the graphicallasso (Friedman et al., 2008), enables convex learning of high-dimensional undirected graphical models.
These methods estimate the inverse covariance of azero-mean Gaussian distribution while penalizing the`1 norm of the o↵-diagonal entries; since the entries
in the inverse covariance correspond to edges in a
Gaussian Markov random field, this method learns asparsely connected graphical model. In recent years,
many algorithms have been proposed for this problem,
Proceedings of the 30 th International Conference on Machine Learning, Atlanta, Georgia, USA, 2013. JMLR:
W&CP volume 28. Copyright 2013 by the author(s).

mwytock@cs.cmu.edu
zkolter@cs.cmu.edu

including projected gradient methods (Duchi et al.,2008), smoothed optimization (Lu, 2009), alternating
linearization methods (Scheinberg et al., 2010), andquadratic approximation (Hsieh et al., 2011).
However, in many prediction tasks we may not want tomodel correlations between input variables. This is thefamiliar generative/discriminative contrast in machinelearning (Ng & Jordan, 2002), where it has been repeatedly observed that in terms of predictive accuracy,
discriminative approaches can be superior (Sutton &McCallum, 2012). This has lead several researchers
within the past year to (independently) propose a generalization of the Gaussian MRF, which we refer to as
the sparse Gaussian conditional random field (CRF),
that allows for discriminative modeling between input
and output variables (Sohn & Kim, 2012), (Yuan &Zhang, 2012), and our own work in (Wytock & Kolter,2012).1 Although previous papers all showed significant promise to the model, they employed o↵-the-shelf
optimization methods (significantly limiting the size of
potential applications) and/or had theoretical resultsthat did not fully highlight the advantages of sparsity.
In this paper, we make three contributions. First, wedevelop a specialized second-order active set method
for estimating sparse Gaussian CRF parameters,
which we show to be several orders of magnitude fasterthan previously proposed algorithms. Second, we develop convergence bounds for the algorithm that establish conditions for exact recovery of underlying models, with rates that specifically highlight the graph degree, improving upon the results in (Yuan & Zhang,2012) in many settings. Third, we present extensive
experimental results on large-scale synthetic data and1

While these formulations were developed independently, they are mathematically identical, and so the model
should be credited to (Sohn & Kim, 2012) as the firstsource, which termed the model “Sparse Conditional Gaussian Graphical Model”. However, in this paper we refer tothe model as the sparse Gaussian CRF, as the discriminative setting coincides with the standard notion of a conditional random field.

Sparse Gaussian conditional random fields

y1

x1

y2

x2

y3

x3

···
···

np + p(p + 1)/2, and thus model can overfit when thenumber of examples m is relatively small.

yp

xn

Figure 1. Illustration of sparse Gaussian CRF model.

two real-world energy forecasting tasks. Here we show
improvement over state-of-the-art methods for wind
power and electrical demand forecasting; these problems are of substantial practical interest, as even small
advances in forecasting accuracy can have notable impact on the efficiency and costs of large power systems.

2. The sparse Gaussian CRF modelLet x 2 Rn and y 2 Rp denote input and output variables for a prediction task. A Gaussian CRF is a loglinear model withp(y|x; ⇤, ⇥) =

1
exp
Z(x)

y T ⇤y2xT ⇥y(1)where the quadratic term models the conditional dependencies of y and the linear term models the dependence of y on x. The model is parameterized by⇤ 2 Rp⇥p , which corresponds to the inverse covariancematrix, and ⇥ 2 Rn⇥p , which maps the inputs to the
outputs; an illustration of the model is shown in Figure 1. Since the CRF is a Gaussian distribution withmean ⇤ 1 ⇥T x, the partition function is given by1= c|⇤| exp
Z(x)xT ⇥⇤ 1 ⇥T x .

(2)

For m data samples, arranged as the rows of X 2
Rm⇥n and Y 2 Rm⇥p , the negative log-likelihoodf (⇤, ⇥) = log p(Y |X; ⇤; ⇥) is given bylog |⇤| + tr Syy ⇤+2Syx ⇥+⇤ 1 ⇥T Sxx ⇥(3)(omitting the constant term c term), where the S
terms are empirical covariancesf (⇤, ⇥) =

Syy =

1 T
1
1
Y Y, Syx = Y T X, Sxx = X T X. (4)mmm

Without regularization, it is straightforward to verify
that this optimization problem is simply a reparameterization of the least squares problem. We can additionally add `2 regularization by adding 2 to the
diagonal elements of S (formally, this corresponds toa Normal-Wishart prior on ⇤ and the columns of ⇥),but again this just corresponds to the regularized leastsquares solution. However, the total number of parameters in this problem (for estimating both ⇥ and ⇤) is

To address this concern, we regularize the maximum
likelihood estimate by adding `1 regularization to ⇥and the o↵-diagonal elements of ⇤; since the `1 norm
encourages sparsity of the parameters, this directly
corresponds to learning a sparse set of edges in ourgraphical model. Our final optimization problem isthen given by minimizing the composite objective
minimize f (⇤, ⇥) + (k⇤k1,? + k⇥k1 )
⇤,⇥

(5)

where k·k1 denotes the elementwise `1 norm, k·k1,? denotes the elementwise `1 norm on o↵-diagonal entries,
and 2 R+ is the regularization parameter.2 This isa convex objective, following from the convexity of the`1 norm and the fact that the log-partition functionof an exponential family graphical model is concave.
Furthermore, the gradients of f are given byr⇤ f (⇤, ⇥) = Syy⇤ 1⇤ 1 ⇥T Sxx ⇥⇤ 1r⇥ f (⇤, ⇥) = 2Sxy + 2Sxx ⇥⇤ 1 ,(6)which in previous work has motivated the use of firstorder non-smooth optimization methods.

3. Optimization
Previous work on the sparse Gaussian CRF (SGCRF)
model has proposed using o↵-the-shelf algorithms tosolve the above optimization problem, including orthantwise quasi-Newton methods (Sohn & Kim, 2012)(specifically, the OWL-QN method of (Andrew & Gao,2007)), and accelerated proximal gradient methods(Yuan & Zhang, 2012) (specifically, the FISTA algorithm of (Beck & Teboulle, 2009)). These methods are
attractive due to their simplicity, and since the gradients can be efficiently computed using (6). Unfortunately, the algorithms still su↵er from relatively slow
convergence (even though they are faster than many
alternative non-smooth first-order methods), and thus
quickly become computationally impractical for large
output and input dimensions.
In this section, we propose a new second-order active
set method for solving the sparse Gaussian CRF. Suchalgorithms have previously been applied to the Gaussian MRF (Hsieh et al., 2011; Olsen et al., 2012), anda general analysis of such methods (showing quadratic
convergence) is presented in (Tseng & Yun, 2009). The
2
It is also possible to introduce di↵erent regularization
parameters for ⇤ and ⇥, though we have found throughour experiments that the optimal settings for these regularization parameters are typically quite similar, so we use
only one for simplicity.

Sparse Gaussian conditional random fields

method here largely mirrors in the approach in (Hsieh
et al., 2011) for the Gaussian MRF, but the precise
formulation is significantly more involved, owing to
the complexity of gradient term of the ⇤ 1 ⇥T Sxx ⇥term in the likelihood. Despite being a second-ordermethod, we show that the resulting algorithm is faster(to reach any accuracy) than previously proposed approaches, and several orders of magnitude faster atachieving solutions to high accuracy.

D ⇤ , D⇥3.1. A second-order active set approach
The basic idea of our method is to iteratively form asecond-order approximation to the objective function(without the `1 regularization term), and then solve an`1 regularized quadratic program to find a regularized
analog of the Newton step. In general notation, tominimize some objective f (x) + kxk1 , we form the
Taylor expansion
f (x+ ) ⇡ g( ) ⌘ f (x)+rx f (x)T

+

1
2

T

r2x f (x)(7)
where rx f (x) and r2x f (x) denote the gradient andHessian respectively. To compute the regularized Newton step d, we solved = arg min g( ) + kx +

k1

(8)

and update the parameters x
x + ↵d for some stepsize ↵, determined by backtracking line search.
In our setting, precise formulations of the gradient andHessian terms are cumbersome, due to the fact that allparameters involved are matrices, but we can concisely
express this second order Taylor expansion using differentials. In particular, (see Appendix A for a fullderivation) the second order Taylor expansion isf (⇤ +⇤, ⇥ +

trSyy
2tr⇤
tr⇤ 1
tr⇤ 1

⇥ ) ⇡ g(⇤ + 2trSyx
1

T

⇥

⇤,
1

tr⇤

⇥ ) ⌘ f (⇤, ⇥) +⇤ +⇥ Sxx ⇥⇤ 1 ⇤ +1
1 T⇥ Sxx ⇥⇤ 1 ⇤ + tr⇤ 1 ⇤ ⇤ 1
⇤⇤
2
T
1
1 T
S
2tr⇤
⇤
⇥ Sxx ⇥ .
xx
⇥
⇤
⇥

⇥ Sxx

⇥

tr⇤

1

T

⇤ +

(9)
As above, we compute the Newton steps D⇤ , D⇥ by
D⇤ , D⇥ = arg min g(
⇤,

⇤,

⇥) +

(10)

⇥

(||⇤ +⇤ k1,? + k⇥ +Algorithm 1 Newton Coordinate Descent for SGCRF
Input: Input features X 2 Rm⇥n and outputs Y 2
Rm⇥p ; regularization parameter
Output: Optimized parameters ⇤, ⇥Initialize: ⇤I, ⇥0
while (not converged) do
1. Determine active sets S⇤ , S⇥ using (11).
2. Find Newton update D⇤ , D⇥ by solving thefollowing optimization using coordinate descent⇥ k1 )

where we use a coordinate descent algorithm to optimize this `1 regularized QP. Since `1 regularization on
the Newton direction tends to push the Newton updates in a direction that increases sparsity and since

arg

⇤,

min
⇥,

(||⇤ +

S̄ =0

g(

⇤,

⇥) +

⇤ k1,? + k⇥ +⇥ k1 ) .

3. Compute a step size ↵ using backtracking line
search, and update⇤⇤ + ↵D⇤ , ⇥⇥ + ↵D⇥ .
end whilethe line search provably converges to step sizes with
↵ = 1 (Tseng & Yun, 2009), the number of nonzero
elements tends to increase as the optimization progresses. For the line search, we additionally need toensure that ⇤ is positive definite, which we ensure by
a common technique of simply defining log |X| to be
infinite if X 6 0. A generic pseudo-code description of
the algorithm is shown in Algorithm 3.1, and a moredetailed presentation is given in Appendix B. Furthermore, a C++ and MATLAB implementation is available at http://www.cs.cmu.edu/~mwytock/gcrf/.
3.2. Computational speedups
In order to make the Newton method efficient numerically, there are a number of needed optimizations.
Again, these mirror similar optimization presented in(Hsieh et al., 2011), but require adaptations for the
CRF case. In practice, the majority of the computational work of the Newton CD method comes fromcomputing the regularized Newton step via coordinate
descent; even though coordinate descent is known tobe an efficient method for solving `1 regularized problems, in our setting we have a total of p(p + 1)/2 + np
di↵erent variables, and it would be infeasible to optimize over them all. Thus, at each iteration of thealgorithm we use an active set method, and only optimize over a variable ( ⇤ )ij or respectively ( ⇥ )ij if|(r⇤ f (⇤, ⇥))i,j | >|(r⇥ f (⇤, ⇥))i,j | >or ⇤ij 6= 0or ⇥ij 6= 0,(11)i.e., if the optimally conditions for that parameter are
violated for the current iterate of the parameters. Because the sparsity resulting from the `1 constraint results in a relatively small active set, this provides aSparse Gaussian conditional random fields

substantial speedup, especially when the optimal solution has high sparsity. We also keep the active set
small by using warm starts; solving the optimization
problem for a decreasing a sequence of the regularization parameter and initializing each successive problem with the previous optimal solution.
Second, in the coordinate descent loop, it is important to cache and incrementally update certain matrix
products, such that we can evaluate subsequent coordinate updates efficiently. This requires that we maintain an explicit form of the matrix products ⇤ ⇤ 1
and ⇥ ⇤ 1 ; crucially, when we update a single coordinate of the ⇥ or ⇤ , we only need to update asingle row of these matrix products, and we can subsequently use only certain elements of these products
to compute each coordinate descent step. Details are
given in Appendix B.
Third, since each step of our Newton method involves
solving an `1 regularized problem itself, it is important that we solve this regularized Newton step only
to an accuracy that is warranted by the overall accuracy of algorithm as a whole. Although more involved
approaches are possible, we simply require that the
inner loop makes no more than O(t) passes over thedata, where t is the iteration number, a heuristic thatis simple to implement and works well in practice.
Last, in cases where nm (which is a setting that
we are crucially interested in for motivating `1 regularization), by not forming the Sxx 2 Rn⇥n matrix
explicitly, we can reduce the computation for products
involving X T X from O(n2 ) to O(mn). Note that thesame considerations do not apply to Syy , since we need
to form an invert the p⇥p matrix ⇤ to compute the gradients. Thus, the algorithm still has complexity O(p3 ),as in the MRF case. However, this highlights anotheradvantage of the CRF over the MRF: when n is large,
just forming a generative model over x and y jointly
is prohibitively expensive. Thus, the sparse Gaussian
CRF significantly improves both the performance andthe computational complexity of the generative model.

4. Theoretical results
As for `1 regularized linear regression and the sparse
Gaussian MRF, it is of significant interest to know
when, if data is generated from a sparse underlying model, the sparse Gaussian CRF is able to recover this model with high probability. In this section, we develop theoretical results that show thesample complexity of the SGCRF grows slower than⌦(d4 (log p + log n)) where d is the maximum degree ofthe output variables in the underlying graph; importantly, this term grows logarithmically in the input andoutput dimensions p and n; relative to the best known
bounds for the special cases of the Gaussian MRF andlinear regression, our bound has a worse dependence
on d, which arises in bounding the error of the Taylor
expansion. This element can likely be improved withmore refined analysis, but we focus here on obtaining
a bound that is logarithmic in p and n, and otherwise
does not depend on on the total number of nonzeros
in the true parameters.
The proof proceeds in the primal-dual witness (PDW)
framework of Wainwright (2009) (that is, we are concerned with the setting of recovering the sparsity ofthe true underlying model) and our analysis mirrors
much of the Gaussian MRF case (Ravikumar et al.,2011); however, as with the optimization, the additional terms in the gradient of the CRF introduce substantial added complexity, which for instance result in
the worse dependence on d. We operate under thefollowing assumptions, similar to assumptions for theGaussian MRF and least-squares settings.
1. True underlying model. The data is generated
according to the modely|x ⇠ N ( ⇤? 1 ⇥? T x, ⇤? 1 )(12)

where each row of [⇤? ⇥? T ] has at most d nonzero
entries (i.e., the vertices corresponding to output
variables in the graphical model of the CRF havemaximum degree d). It is straightforward to generalize this to the case of sub-Gaussian noise, but
we assume the Gaussian model for simplicity.
2. Column normalization. The columns of the
input feature matrix havep bounded `2 norm suchthat maxj=1,...,n kXj k2 / m  cX . This same assumption is used in the corresponding analysis ofthe `1 regularized least-squares case. Importantly,
in the CRF case we make no assumptions about
the distribution of x.
3. Restricted convexity. Letting Si denote thenonzero indices of the ith column of ⇥? (i.e., theedges between inputs and the ith output), we have
✓
◆
1 T
X
X
>0
(13)
min
S
m Si ii.e., this term is strictly convex when restricted
to the true active set. This is again a common
assumption for `1 methods, but note that we do
not require a restricted convexity assumption on⇤? as the logdet term is already strictly convex.

4. Mutual incoherence. This is the most subtleof the assumptions, and one which can often be
violated in practice, yet it is required for the PDW

Sparse Gaussian conditional random fields

approach. Denoting the Hessian of the objective
as H = r2⇤,⇥ f (⇤, ⇥) and S the set of all nonzero
entries of ⇤? and ⇥? , we require that|||HS̄S (HSS ) 1 |||1  1

↵

(14)

for some ↵ > 0 where where |||·|||1 denotes thematrix infinity norm, the maximum absolute row
sum. This condition stipulates that the edges inthe true active set are not too correlated withedges outside, and mirrors the same assumption for the least-squares and Gaussian MRF approaches (though of course with di↵erences owing
to the precise form of the Hessian). We give an
illustrative example of this condition for simple
graphs in Appendix D.
Theorem 1. Under the above assumptions, given asample size and regularization parametermc1 d4 (1 + 8/↵)2 log(pn)rlog(pn) + log 4
1
c2 ↵
m

then with probability at least 1(15)

c3 exp( c4 m 2 )1. The solution to the `1 regularized optimization˜ ⇥,˜ have nonzero entries that are a
problem, ⇤,strict subset of the nonzero entries of ⇤? , ⇥? .
2. The solution satisfies the elementwise bounds
˜
max{k⇤˜ ⇥? k1 } ⇤? k1 , k⇥rlog(pn) + log 4
c5 (1 + 8↵ 1 )m(16)

where c1 , . . . , c5 denote constant terms.
The proof of this theorem is quite lengthy and deferredto Appendix C (where we also provide an explicit definition of the constants and a lengthier definition of theassumptions). Intuitively, this theorem shows that asample size of m = ⌦(d4 (log p + log n)) is sufficientto guarantee with high probability that solving the `1
regularized MLE recovers a subset of the true edge
structure, and that the recovered parameters are close
to the true parameters. Note that in many settings
this is an improvement over the bound in (Yuan &Zhang, 2012), which e↵ectively requires a sample size⌦(s(log p+log n)) where s is the total number of edges
in the graph; for graphs with a fixed low degree (suchas a chain grain) s can grow linearly in p or n, whereasd remains constant, and so this represent a significant
improvement—indeed, as we show in our experimental
results, the empirical scaling does indeed depend only
logarithmically on p, even if s increases linearly in p.3
3

The bounds are not directly comparable, since Yuan

5. Experimental results
Here we experimentally evaluate several aspects of the
proposed model and algorithm on synthetic data andtwo real-world energy forecasting problems, the tasksof predicting upcoming wind power from multiple windfarms and the task of predicting upcoming electricaldemand over multiple utility zones. For the latter twocases, we demonstrate state-of-the-art results. The
wind prediction task is from the 2012 Global Energy
Forecasting Competition (Hong, 2012), a contest recently held on Kaggle to forecast wind power; here
our algorithm improves upon our own submission tothis contest by 5.5% (our entry was a top-5 entry
that used least-squares with the same features and was
2.5% worse than the winning entry). For load forecasting, we use real-world load data from the PJM system
operator (available at http://www.pjm.com/) and improve upon the deployed PJM forecasts by 19%. We
also highlight the performance of the algorithm relative to its theoretical bounds, and the optimization
performance of our Newton method (which in all cases
substantially improves upon previous methods).
5.1. Synthetic data
Exact subset recovery. Our first experiments illustrate when the model is able to exactly recover theunderlying graph structure of a true model, and illustrates that the overall dependence given in our theory
looks to fit the observed results. Specifically, we generate data from a chain CRF, where each output variable
is connected to two others, and each input variable isconnected to one output. To represent the chain we
use the true parameter ⇤? with ⇤?ii = 1 on the diagonal, ⇤?ij = 0.2 on the super diagonal and a diagonal⇥? with ⇥?ii = 0.2.
In Figure 2 (top), we vary m for di↵erent choices of pand observe that once m passes a certain threshold werecover the support of the true parameters with highprobability—scaling the x-axis by log p demonstrates
the same theoretical dependence on p as shown in ourtheory. Importantly, note that in this case, the total number of edges in the graph, s, increases linearlyin p whereas the maximum vertex degree is fixed at
d = 3. Thus, our bound captures the overall scalingof the model, whereas the bound of (Yuan & Zhang,2012) would be significantly looser in this case. For the
Figure 2 (bottom), we increase n by adding irrelevant
features (features that are not connected to the output& Zhang (2012) bounds only the Frobenius norm, and requires a looser restricted isometry property. Nonetheless,we see no direct way of providing a dependence on graphdegree using the analysis methods in this past work, so thisrepresents a notable improvement.

4

1

10

0.8

10

0.6

10

0

0.4

p=50
p=100
p=200
p=400

0.2

0

1000 2000 3000 4000 5000
m

500
1000
m / log p

−2

10

−4

10

−6

10

0

−8

10

1500

0

1000 2000 3000 4000 5000 6000
Time (seconds)

Figure 3. Suboptimality of solution versus time for Newton CD versus previously considered algorithms for sparse
Gaussian CRF, OWL-QN and FISTA.

1
0.8

1

0.6

SGCRF
SGMRF

0.9

0.4

n=50
n=100
n=200
n=400

0.2

0.8
MSE

Probability of success

Newton CD
OWL−QN
FISTA

2f − f*

Probability of successSparse Gaussian conditional random fields

0
0

1000

2000
m

3000

200
400
m / log n

600

Figure 2. (top) Fraction of 20 trials in which support ofthe estimated parameters match that of the true parameters, increasing n, p; (bottom) adding irrelevant features,increasing n.

variables); again, we observe a logarithmic dependence
on the input dimension.
Optimization performance. Because the chain
CRF is a rather limited example, for the remaining
synthetic examples we generate data from more complex model. In particular, follow a similar procedure
as in (Yuan & Zhang, 2012), and generate ⇤ and ⇥with 5(n + p) random unity entries (the rest beingzero), and set the diagonal of ⇤ such that the condition number of ⇤ equals n + p. We sample x from azero-mean Gaussian with full covariance, square half
the entries, and then normalize the columns to haveunit variance. We use this same process for the next
three experiments, but vary problem size to make theexperiments computationally feasible in all cases.
Figure 3 shows the suboptimality of each method interms of the objective function f f ? (where f ? is computed by running our Newton CD approach to numerical precision) versus execution time on a 2.4GHz Xeon
processor; this problem has size p = 1000, n = 4000,
and m = 2500. On this problem the Newton CD approach converges to high numerical precision within
about 81 seconds, while FISTA and OWL-QN still
don’t approach this level of precision after two hours.
It is also important to note that the Newton CD approach also reaches all intermediate levels of accuracy
faster than the alternative approaches, so that the al-

0.7
0.6
0.5
0.4
1
10

0

10

−1

10
λ

−2

10

−3

10

Figure 4. Generalization performance (measured by meansquared error of the predictions) for the Gaussian MRF
versus CRF, with problem size n = 200, p = 50, m = 50.

gorithm is preferable even if only intermediate precision is desired. Indeed, we note previous works (Sohn& Kim, 2012; Yuan & Zhang, 2012) considered maximum problem sizes of np ⇡ 105 due to the time required for training; since much of the appeal of `1 approaches lies precisely in the ability to use large featuresizes, this has significantly limited the applicability of
the approach. We thus believe that our proposed algorithms opens the possibility of substantial new applications of this sparse Gaussian CRF model.
Comparison to MRF. Our next experiment compares the discriminative CRF modeling to a generative
MRF model. In particular, an alternative approach to
our framework is to use a sparse Gaussian MRF to
jointly model x, y as a Gaussian, then compute y|x.
Figure 4 shows the performance of the Gaussian MRF
versus CRF, measured by mean squared error in thepredictions on a test set, over a variety of di↵erent
parameters. The CRF substantially outperforms
the MRF in this case, due to two main factors: 1)
the x variables as generated by the above process are
not Gaussian, and thus any Gaussian distribution willmodel them poorly; and 2) the x variables are correlated and have dense inverse covariance, making it
difficult for the MRF to find a sparse solution.

Sparse Gaussian conditional random fields
0.8
SGCRF
LS

MSE

0.6
0.4
0.2
0

100

200
300
400
Training examples (m)500

Figure 5. Generalization performance (MSE for the bestchosen via cross-validation), for the sparse Gaussian CRF
versus `2 regularized least squares. Here n = 800, p = 200.

Finally, we note again that in addition to the performance benefits, the CRF has substantial computational benefits. Modeling x and y jointly requires
computing and inverting their joint covariance, whichtakes time O((n + p)3 ); in contrast, the corresponds
operations for the CRF case are O(np3 ), which is substantially faster for even modestly large n. Indeed, forthe two real-world experiments below, we were unable
to successfully optimize a joint MRF using the QUIC
algorithm of (Hsieh et al., 2011) (itself amongst thefastest for solving the sparse Gaussian MRF), after
running the algorithm for 20 hours.
Sample size. Finally, to illustrate the benefit of `1
regularization over traditional (`2 regularized) multiple least-squares estimation, we evaluate generalization performance versus sample size, shown in Figure 5. This figure shows performance measured bymean squared error of the `1 regularized sparse Gaussian CRF versus traditional least-squares with `2 regularization; here, for each m we choose the `1 and `2
regularization parameters using a cross validation set,
then evaluate the MSE on a test set. As the sample
size increases, the performance of the two methods becomes similar (in the limit of infinite data with fixedn and p, they will of course be equivalent); however,as expected, for small samples sizes the `1 regularization method performs much better, being able to takeadvantage of the sparsity in the underlying model.
5.2. Application to energy forecastingWind power forecasting. We here apply the sparse
Gaussian CRF model to the wind power forecastingtask from the Global Energy Forecasting 2012 competition (Hong, 2012), a recent Kaggle competition
for predicting upcoming wind power at seven di↵erent nearby wind farms for a time horizon of 48 hours.
The input data for this problem consisted of previous
power outputs for the wind farms (going as far backas the past 36 hours), and wind speed forecasts for

Figure 7. Sparsity patterns of estimated ⇤ and ⇥ parameters for the wind forecasting task. White denotes zero
values, and a wind farms are grouped together in blocks.

the upcoming 48 hours. From this input we generated features that consisted of: 1) the past 8 hoursof power for each wind farm, and 2) 10 RBF features
placed around each forecasted wind speed, to capture
non-linear dependencies on the wind speed itself. In
total, this lead to p = 336 dimensional outputs andn = 3417 dimensional inputs. We heavily optimized
these features for the competition, and using these features with ordinary least-squares resulted in a top-5
finish in the competition (out of 134 entrants).
Figure 6 shows the performance of the sparse Gaussian
CRF on the wind forecasting task, analyzed across several dimensions. First, the figure on the left shows performance of method for varying ; also shown in thebest performance of `2 regularized least-squares. For
properly chosen , the algorithm outperforms leastsquares (using the exact same features), by 5.5%. Fora domain such as wind power forecasting, where there
is a limit to the possible performance (wind is an inherently stochastic phenomenon, so exact forecasts are
not possible), and since the least-squares solution inthis case is already using highly optimized features,
this represents a substantial improvement. The di↵erence in performance become even more pronounced forsmaller sample sizes, as shown in the Figure 6 (center),which shows MSE (using chosen by hold out crossvalidation), for a variety of sample sizes. Finally, tohighlight the optimization performance on real data
as well, we shown in Figure 6 (right) the optimization
objective versus training time for the di↵erent optimization algorithm. Again, the Newton CD algorithm
vastly outperforms FISTA and OWL-QN, converging
to high accuracy after 160 minutes, whereas the latterdo not reach reasonable accuracy after several hours.
Finally, a significant advantage of the sparse Gaussian CRF approach is that the sparsity pattern of theresulting model can be interpreted directly as conditional dependencies between variables, and thus thesparsity pattern itself can be very informative. Figure
7 shows the sparsity patterns in ⇤ and ⇥ for the windforecasting task; they illustrate a clear temporal andspatial dependence between the di↵erent wind farms.
Electrical demand forecasting. We further applySparse Gaussian conditional random fields
4

0.4
SGCRF
LS

0.5

10
SGCRF
LS

0.38

2

10

0

0.4

10
f − f*

MSE

MSE

0.36
0.34

−4

0.32

0.3

−2

10
10

Newton CD
OWL−QN
FISTA

−6

0.3

10

−8

0.2
0
10

−1

−2

10
λ

10

200
400
600
800
Training examples (m)

10

0

100
200
300
Time (minutes)

400

Figure 6. Performance of SGCRF on wind power forecasting showing (left) generalization performance for varying valuesof with vertical line denoting value chosen by cross-validation, (center) performance versus least-squares on di↵erent
sample sizes, and (right) optimization performance of the Newton CD approach versus alternatives.

shows that by jointly predicting over all the zones, weare able to improve substantially upon PJM’s alreadystate-of-the-art point forecasts.

0.12
SGCRF
PJM forecast

MSE

0.1

6. Conclusion and discussion

0.08
0.06
0.04

−1

−2

10

10

−3

10

λ

Figure 8. Generalization performance for forecasting future demand for 24 hours, compared under MSE to PJM’s
own forecasts with vertical line denoting chosen by crossvalidation.

the model to the task of predicting future electricaldemand for zones operated by PJM (a system operator
for coordinating electricity generation and delivery forseveral Eastern U.S. states). In particular, the goal is
to forecast upcoming electrical demand for the next 24hours over 15 di↵erent zones in the system.
Electricity forecasting is a well-studied problem (Soliman & Al-Kandari, 2010), and PJM already employs asophisticated forecasting system in its operation (Various, 2012) to predict a subset of the zones; rather thantry to build an entirely new forecast, we use these previous point forecasts as input features (along with pastenergy consumption and time-of-day features) to predict future demand. The goal is thus to use a combination of the existing predictions to predict even moreaccurately, and if we can improve upon the PJM forecasts, this means that we are e↵ectively combining existing information to ultimately deliver a better prediction. For this problem, the dimension of input andoutput are p = 350 and n = 860. We present theseresults here more briefly, but the key performance element we want to highlight is in Figure 8; the figure

The sparse Gaussian conditional random field enjoys
many benefits of existing methods for learning highdimensional Gaussian graphical models; we believethat the advances put forward in this paper make themodel significantly more practical for large-scale problems, and also significantly advance our theoretical understanding of the method. Furthermore, the empirical results presented here on wind power and demand
forecasting are of substantial practical interest, andthe improvements presented here have the potential
for notable e↵ects on power system efficiency.
Two future directions seem particularly promising.
First, it would be worthwhile to use regret-based approaches to develop alternate convergence rates under
weaker assumptions than those we use. Although exact feature selection is not possible even for the leastsquares case when inputs are very highly correlated, it
is nonetheless possible to obtain regret bounds thatbound the loss versus that of the true model, e.g.
(Bartlett et al., 2012); such directions are likely to be of
substantial interest here, since we do expect to often bein situations where input features are correlated. Second, from an application standpoint in energy systems
in particular, there are a huge number of forecasting
problems that share similar properties as wind powerand demand; of crucial importance, however, is developing control algorithms that can exploit these probabilistic forecasts in the planning stage. Developing
such algorithms will allow high-dimensional graphicalmodels such as the Gaussian CRF to have an immediate impact on these globally important domains.

Sparse Gaussian conditional random fields

References
Andrew, Galen and Gao, Jianfeng. Scalable training ofl 1-regularized log-linear models. In Proceedings ofthe 24th international conference on Machine learning, pp. 33–40. ACM, 2007.
Banerjee, O., Ghaoui, L. El, and d’Aspremont, A.
Model selection through sparse maximum likelihood
estimation for multivariate gaussian or binary data,
2008.
Bartlett, Peter L, Mendelson, Shahar, and Neeman,
Joseph. `1 -regularized linear regression: Persistence
and oracle inequalities. Probability theory and related fields, pp. 1–32, 2012.
Beck, Amir and Teboulle, Marc. A fast iterative
shrinkage-thresholding algorithm for linear inverse
problems. SIAM Journal on Imaging Sciences, 2(1):183–202, 2009.
Duchi, J. C., Gould, S., and Koller, D. Projected subgradient methods for learning sparse Gaussians. In
Proceedings of the Conference on Uncertainty in Artificial Intelligence, 2008.
Friedman, J., Hastie, T., and Tibshirani, R. Sparse inverse covariance estimation with the graphical lasso.
Biostatistics, 9(3):432–441, 2008.
Hong, T. Global energy forecasting competition, 2012.
URL http://www.gefcom.org.
Hsieh, C.-J., Sustik, M. A., Dhillon, I. S., and Ravikumar, P. Sparse inverse covariace matrix estimation
using quadratic approximation. In Neural Information Processing Systems, 2011.
Lu, Z. Smooth optimization approaches for sparse inverse covariance selection. SIAM Journal on Optimization, 19(4):1807–1827, 2009.
Ng, A. Y. and Jordan, M. I. On discriminative vs. generative classifiers: a comparison of logistic regression
and naive Bayes. 2002.
Olsen, Peder A, Oztoprak, Figen, Nocedal, Jorge, andRennie, Stephen J. Newton-like methods for sparse
inverse covariance estimation. Optimization Online,
2012.
Ravikumar, Pradeep, Wainwright, Martin J, Raskutti,
Garvesh, and Yu, Bin. High-dimensional covariance estimation by minimizing 1-penalized logdeterminant divergence.
Electronic Journal of
Statistics, 5:935–980, 2011.

Scheinberg, K., Ma, S., and Goldfarb, D. Sparse inverse covariance selection via alternating linearization methods. In Neural Information Processing
Systems, 2010.
Sohn, Kyung-Ah and Kim, Seyoung. Joint estimation of structured sparsity and output structure
in multiple-output regression via inverse-covarianceregularization. In Proceedings of the Conference on
Artificial Intelligence and Statistics, 2012.
Soliman, S. A. and Al-Kandari, A. M. Electrical Load
Forecasting: Modeling and Model Construction. Elsevier, 2010.
Sutton, C. and McCallum, A. An introduction to conditional random fields. Foundations and Trends in
Machine Learning, 4(4):267–373, 2012.
Tseng, Paul and Yun, Sangwoon. A coordinate gradient descent method for nonsmooth separable minimization. Mathematical Programming, 117(1):387–
423, 2009.
Various. PJM Manual 19: Load Forecasting and Analysis. PJM, 2012. Available at: http://www.pjm.
com/planning/resource-adequacy-planning/~/
media/documents/manuals/m19.ashx.
Wainwright, M.J.
Sharp thresholds for highdimensional and noisy sparsity recovery using `1 constrained quadratic programming (Lasso). IEEE
Transactions on Information Theory, 55(5):2183–
2202, 2009.
Wytock, Matt and Kolter, J. Zico. Sparse conditional
gaussian random fields. In NIPS Workshop on Log
Linear Models, 2012.
Yuan, Xiao-Tong and Zhang, Tong. Partial gaussian
graphical model estimation. CoRR, abs/1209.6419,
2012.

Appendices to “Sparse Gaussian Conditional Random Fields:
Algorithms, Theory, and Application to Energy Forecasting”

A. Derivation of Gradient, Hessian, and Di↵erentials
Here we derive analytic expressions for the gradient, Hessian, and various di↵erentials of the log likelihood
function. Recall that the log-likelihood is given byf (⇤, ⇥) =log |⇤| + tr⇤Syy + 2tr⇥T Sxy + tr⇤ 1 ⇥T Sxx ⇥(1)

We adopt the di↵erential matrix calculus notation from (Magnus & Neudecker, 1988) where for a matrix A 2
Rm⇥n and a function f : Rm⇥n ! R, dk f (A; U ) : Rm⇥n ⇥ Rm⇥n ! R denotes the k-th di↵erential of the function
f evaluated at U . For example, the first di↵erential can easily be expressed in terms of the gradientdf (A; U ) = trrA f (A)T U.

(2)Other derivatives for functions of matrices (i.e., Hessians or higher order terms) are cumbersome to represent
directly, but the di↵erentials themselves can typically be expressed compactly; indeed, it is often simplest to
first derive these di↵erentials and then use them to determine analytical expressions for the Hessians and higher
order derivatives. Furthermore, the Taylor expansion of a function can be represented directly in terms of its
di↵erentials; for instance the second order approximation is given byf (A +) ⇡ f (A) + df (A;1) + d2 f (A;2) ⌘ f (A) + vec(rA f (A))T vec( ) +1vec( )T (r2A f (A)) vec( ) (3)2where vec denotes the vectorization of a matrix (forming a vector by concatenating the columns), and r2A f (A)denotes the Hessian.
Using standard rules of di↵erential calculus, we can compute the first and second order di↵erentials of the
log-likelihood f (⇤, ⇥),df (⇤, ⇥; U, V ) = trSyy U + 2trSyx Vtr⇤ 1 U + 2tr⇤ 1 ⇥T Sxx V

tr⇤ 1 ⇥T Sxx ⇥⇤ 1 U(4)

and from this expression we can easily determine the relevant gradientsr⇤ f (⇤, ⇥) = Syy⇤ 1⇤ 1 ⇥T Sxx ⇥⇤ 1

(5)r⇥ f (⇤, ⇥) = 2Sxy + 2Sxx ⇥⇤ 1 .
Similarly, we can di↵erentiate again to find the second di↵erential
d2 f (⇤, ⇥; U, V ) = 2tr⇤ 1 U ⇤ 1 ⇥T Sxx ⇥⇤ 1 U + tr⇤ 1 U ⇤ 1 U + 2tr⇤ 1 V T Sxx V

4tr⇤ 1 U ⇤ 1 ⇥T Sxx V. (6)Combining the first and second di↵erential gives the full second order Taylor expansion shown in the paper. It
also lets us determine the Hessian itself, which we use in the incoherence condition for the theoretical results
 1⇤ ⌦ (⇤ 1 + 2⇤ 1 ⇥T Sxx ⇥⇤ 1 )2⇤ 1 ⌦ ⇤ 1 ⇥T Sxx
r2⇤,⇥ f (⇤, ⇥) =(7)1
1
2⇤ ⌦ Sxx ⇥⇤2⇤ 1 ⌦ Sxx

B. Detailed Description of Newton Coordinate Descent Method
We present a detailed description and full pseudo-code for the Newton coordinate descent algorithm. The
derivation mirrors that in (Hsieh et al., 2011). The complete method is shown in Algorithm 1, with the coordinate

Appendices to “Sparse Gaussian conditional random fields”

Algorithm 1 Newton Coordinate Descent for SGCRF
Input: Input features X 2 Rm⇥n and outputs Y 2 Rm⇥p ; regularization parameter ; algorithm parameters✏, , ↵, .
Output: Optimized parameters ⇤, ⇥Initialize: ⇤
I, ⇥
0, ⌃
⇤ 1
while (not converged) do
1. Compute the gradient, determine active sets S⇤ , S⇥ using (14), and check for convergence.
2. Find regularized Newton direction D⇤ , D⇥ using Algorithm 2.
3. Initialize ↵1 and compute(trr⇤ f (⇤, ⇥)Tµ⇤ + trr⇥ f (⇤, ⇥)T⇥ + k⇤ +⇤ k1? + k⇥ +⇥ k1 ).

while (insufficient descent) do
1. Compute the Cholesky decomposition LLT = ⇤ + ↵D⇤ , continuing if not positive definite
2. Check descent f (⇤ + ↵D⇤ , ⇥ + ↵D⇥ ) < f (⇤, ⇥) + ↵ µ and break if satisfied
3. ↵↵end while
end while
descent inner loop for computing the approximation to the Newton direction given in Algorithm 2. This processrepeats until the solution converges to a within a specified tolerance, checked using the KKT conditions.
Next, we derive the coordinatewise updates for the inner loop and highlight the key optimizations that are usedin order to achieve fast performance.
B.1. Coordinate descent updates for the Newton approximation
To begin, note that for a fixed ⇤ and ⇥, the regularized Newton direction is given by the solution to thesecond-order Taylor expansion which for our problem has the form
h(

⇤,

⇥) =

trSyy

⇤ + 2trSyx

tr⇤ 1 ⇥T Sxx ⇥⇤ 1(||⇤ +⇥tr⇤ 1⇤ + tr⇤⇤ k1,? + k⇥ +1⇤

+ 2tr⇤ 1 ⇥T Sxx⇤⇤

1⇥T Sxx ⇥⇤ 1

⇥
⇤

1
tr⇤ 1
2

⇤⇤

1

⇤+

(8)

⇥ k1 )

We split the updates into three cases. First, we consider optimizing over a diagonal element of D⇤ by findingµ = arg minµ h( ⇤ + µei eTi , ⇥ ) which has the explicit form
minimizeµ1 2⇥ 2µ ⌃ii + 2⌃ii
2|⇤ii + µ|where ⌃ = ⇤ 1 and

ii

⇤

+µ

⇥

⌃ii + (Syy )ii

ii + (⌃U ⌃)ii⇤
2(⌃⇥T Sxx V ⌃)ii + 2( U ⌃)ii +(9)= ⌃⇥T Sxx ⇥⌃.

Next, note that for two symmetric matrices A, B, not necessarily equal, the symmetric update is given byarg min trA(U + µ(ei eTj + ej eTi ))B(U + µ(ei eTj + ej eTi ))µ= arg min µ2 trA(ei eTj + ej eTi )B(ei eTj + ej eTi ) + µtrAU B(ei eTj + ej eTi ) + µtrA(ei eTj + ej eTi )BUµ(10)= arg min µ2 (Aii Bjj + 2Aij Bij + Ajj Bii ) + 2µ((AU B)ij + (AU B)ji )µApplying this equivalence twice, once with A = B = ⌃ and again with A = ⌃, B =the the symmetric update

Appendices to “Sparse Gaussian conditional random fields”

Algorithm 2 Coordinate descent inner loop
Input: S empirical covariance, regularization parameter, S⇤ , S⇥ active sets and ⇤, ⇥ current estimates
Output: Approximate regularized Newton direction D⇤ , D⇥
Initialize: D⇤
0, D⇥
0, U
0, V
0
while (not converged) dofor coordinate (i, j) in S⇤ do
1. Find µ by solving (9) or (11), using U and V for efficiency.
2. Symmetrically update D⇤ and two rows of U(D⇤ )ij , (D⇤ )ij(D⇤ )ij + µ

Ui

Ui + µ⌃j

Uj

Uj + µ⌃i

where ⌃i denotes the ith row of ⇤ 1 .
end forfor coordinate (i, j) in S⇥ do
1. Find µ by solving (12), using U and V for efficiency.
2. Update D⇥ and one row of V
(D⇥ )ijVi(D⇥ )ij + µ
Vi + µ⌃j

end forend whileTT⇤ + µ(ei ej + ej ei ),for an o↵-diagonal element of matrix D⇤ , µ = arg minµ h(⇥
minimize µ2 ⌃2ij + ⌃ii ⌃jj + ⌃iiµ

2µ [ ⌃ij + (Syy )ij

jj + 2⌃ij

ij + ⌃jj

ij + (⌃U ⌃)ij

ij

ii

⇤

⇥ ) is given by+

ji + (U ⌃)ij + ( U ⌃)ji ] +(11)2 |⇤ij + Uij + µ|where

= ⌃⇥T Sxx V ⌃. Finally, we consider optimizing over an element of D⇥
minimize µ2 [⌃jj (Sxx )ii ] + µ [2(Sxy )ij + 2(Sxx ⇥⌃)ij + 2(Sxx V ⌃)ijµ2(Sxx ⇥⌃U ⌃)ij ] +(12)|⇥ij + Vij + µ|
Each equation has a quadratic form and thus can be solved in closed form. The second two have an `1 penalty
and the form minµ 21 aµ2 + bµ + |c + µ| which has the solution✓◆bµ = c + S /a c(13)aB.2. Optimizations
As in the case of the MRF (Hsieh et al., 2011), there are several modifications to a naive solution that significantly
reduce the running time of the algorithm.
First, consider the matrix products involved in the coordinatewise updates above. A naive implementation of thecoordinate descent algorithm would require O((n + p)2 ) operations even though the majority of the elements areunchanged from one iteration. However, by caching products of the static matrices and maintaining a factorized
form of the products involving ⇤ and ⇥ , specifically U = ⇤ ⌃, V = ⇥ ⌃, we reduce this to O((n + p)). As aconsequence, at each iteration of the loop we must update the rows of U and V corresponding to the coordinatesof ⇤ and ⇥ .

Appendices to “Sparse Gaussian conditional random fields”

Next, we describe how we drastically reduce the coordinate descent active set. At each iteration of the outer loop,
we fix the active set using the current nonzero coordinates and the KKT conditions of the objective function.
We include a coordinate of ⇤, respectively ⇥, if|(r⇤ f (⇤, ⇥))i,j | >|(r⇥ f (⇤, ⇥))i,j | >or ⇤ij 6= 0or ⇥ij 6= 0.

(14)

Since the size of this active set is determined by the number of nonzero elements in the parameters, for sparse
solutions the speed up is very significant. Note that although we fix the active set before beginning coordinate
descent, as in the MRF case (Hsieh et al., 2011), we still have convergence guarantees for the overall algorithm.
Finally, note that when taking a step we must ensure sufficient descent and that the ⇤ parameter remains inthe semidefinite cone. We accomplish this using the Cholesky decomposition, which is also used for efficiently
computing ⇤ 1 .

C. Theoretical Analysis
We will make the following assumptions about the input and output variables X and Y . These mirror similarassumptions in (Wainwright, 2009) and (Ravikumar et al., 2011), and we will discuss the precise di↵erences.
First, the analysis here proceeds on the assumption that there is a true underlying model generating the test
data, of the prescribed form (i.e., the data is generated according to a sparse Gaussian CRF). It is trivial toextend this analysis to the case of sub-Gaussian noise, but we simply assume Gaussian noise for simplicity ofpresentation
Assumption 1. Underlying model The data is generated according toy|x ⇠ N ( ⇤? 1 ⇥? T x, ⇤? 1 ).

(15)

?T

where each row of [⇤? ⇥ ] has at most d nonzero entries (i.e., the vertices corresponding to output variables in
the graphical model of the CRF have maximum degree d).
For simplicity, we will also denote ⌃? = ⇤? 1 .
Our second assumption is a restricted convexity requirement, which ensures that the optimization problem
restricted to the active set is unique. This is a common assumption for `1 approaches (the same condition
appears in the least-squares analysis of (Wainwright, 2009)), and the only extension here is that we require thisto hold for each output variable.
Assumption 2. Restricted convexity For each output i, let Si denote {j : ⇥ji 6= 0} (i.e., Si is the “active
set” of edges directly connecting an input to yi ), we have that
✓
◆
1 T
X
X
> 0.
(16)
min
S
m Si iThe next assumption is more subtle (and quite strict in practice), but is again typical for exact subset selection
proofs for `1 approaches. Namely, we require a mutual incoherence assumption, which e↵ectively ensures thatthe connections in the CRF that correspond to the “true” edges do not correlate too much with edges that are
not the support set.
Assumption 3. Mutual incoherence Let S denote the active set of all variables in vector formvec(supp{⇤? })S=
vec(supp{⇥? })(17)

where supp denotes the support function (the indicator of whether an element is nonzero), and let S̄ denote its
complement. Then for H = r2⇤,⇥ f (⇤, ⇥) defined above|||HS̄S (HSS ) 1 |||1  1↵

for some ↵ > 1, where |||·|||1 denotes the matrix infinity norm, the maximum absolute row sum.

(18)

Appendices to “Sparse Gaussian conditional random fields”

Our first lemma shows that the gradients r⇤ f (⇤? , ⇥? ) and r⇥ f (⇤? , ⇥? ) (the gradients evaluated at the true
parameters) are small (in infinity norm) with high probability given samples on the order of m = ⌦(log n + log p).
The proof (shown below) follows from a standard bound on Gaussian random variables, and from Lemma 1 in(Ravikumar et al., 2011).
Lemma 1. Given data generated by the model in Assumption 1 we have that
⇢
m✏2
?
?
(19)
P (kr⇥ f (⇤ , ⇥ )k1 > ✏)  2np exp
8c2 ? c2X
pwhere c ? = maxi ⌃?ii , and cX = maxj=1,...,n kXj k2 / m; the maximum normalized `2 norm over columns of X.
Furthermore,⇢m✏2
P (kr⇤ f (⇤? , ⇥? )k1 > ✏)  4p2 exp(20)3200c2 ?
for 0 <

< 40c ? .

The next lemma is a generic primal-dual witness approach, mirroring exactly the derivation in (Wainwright,2009), but presented in a generic form. For the presentation here, we will use a generic optimization problem
minimize f (✓) + k✓k1 , though we will apply this specifically to our CRF problem momentarily. Intuitively, thelemma states conditions for which optimizing over the known support set is equivalent to optimizing with the `1
penalty.
Lemma 2. Consider some sparse ✓? with S = supp(✓? ), and consider the two optimization problems✓ˆ = arg min f (✓) + k✓k1✓

(21)

✓˜ = arg min f (✓) + k✓k1 .
✓,✓S̄ =0

Define= ✓˜✓? , and R( ) =˜ + r✓ f (✓? ) + r2 f (✓? ) . Then if the following conditions holdr✓ f (✓)✓

1. The solution ✓˜ is unique.
2. Mutual incoherence holds, i.e., |||(r2✓ f (✓? ))S̄S (r2✓ f (✓? ))SS1 |||1  1↵

3. max{kr✓ f (✓? )k1 , kR( )k1 }  ↵/8ˆthen the `1 solution recovers the restricted solution, ✓˜ = ✓.
In our setting, theand R( ) are themselves matrices and thus slightly more complex. Thus, for subsequentlemmas, we define the following terms that we use to quantity the error of the second order Taylor expansionsof our particular log-likelihood, evaluated at the true parameters. For any ⇤, ⇥, we define ⇤ = ⇤ ⇤? and⇥? and
⇥ =⇥

⇤
⌘
.
(22)
⇥
Define
R⇤ (
R⇥ (

⇤,
⇤,

?
?
⇥ ) = r⇤ f (⇤ , ⇥ )?
?
⇥ ) = r⇥ f (⇤ , ⇥ )r⇤ f (⇤? +?

r⇥ f (⇤ +

⇤, ⇥

?

+

?

+

⇤, ⇥

?
?
⇥ ) + d(r⇤ f (⇤ , ⇥ );?
?
⇥ ) + d(r⇥ f (⇤ , ⇥ );

⇤,

⇥)

⇤,

⇥ ),

(23)

which are the residuals of the first order Taylor expansion of the gradient (i.e., the errors in the second order
Taylor expansion of the function itself), and
R⇤ ( ⇤ , ⇥ )R( ) =.
(24)
R⇥ ( ⇤ , ⇥ )The next lemma bounds the residual kR( )k1 in terms of the distance from the true parameters, k k1 .

Appendices to “Sparse Gaussian conditional random fields”

Lemma 3. Under the definitions above, if
k k1 

1
min
d

⇢

1 c⇥?
,
3c⌃? 2(25)

thenkR( )k1  206c4⌃? c2⇥? c2X d2 k k21(26)

where c⌃? = maxi,j ⌃?ij and c⇥? = maxi,j ⇥?ij .

Finally, we show that when the gradient evaluated at the true model are sufficiently small, then k k1 itself issmall.
Lemma 4. Under the model above, suppose thatmax{kr⇤ f (⇤? , ⇥? )k1 , kr⇥ f (⇤? , ⇥? )k1 } 


⇢
1
1
1
min
,
2cH ?
3c⌃? d 412c4⌃? c2⇥? c2X d2

.

(27)

Then
k k1  2cH ? (max{kr⇤ f (⇤? , ⇥? )k1 , kr⇥ f (⇤? , ⇥? )k1 } + )(28)?
where cH ? = maxi,j Hij, the maximum element of the Hessian evaluated at the true parameters.

These elements allow us to prove the desired theorem.
Theorem 1. Using assumptions 1-3 above, suppose we have sample sizem4122 C 2 d4 (1 + 8/↵)2 log(pn)where C = max{3c⌃? , c⇥⇤1 , c4⌃? c2⌃? c2X } and chooseasp(8/↵)c ? cX 3200
then with probability greater than 1

(29)

r

log(pn) + log 4m(30)

c1 exp( c2 m 2 ) we have˜ ⇥,˜ has nonzero entries that are a strict subset
1. The solution to the `1 regularized optimization problem, ⇤,?
?
of the nonzero entries of ⇤ , ⇥2. The solution satisfies the elementwise bounds
˜
max{k⇤˜⇤ k1 , k⇥?

?

⇥ k1 }  2(1 + 8↵

Proof. Let

1

)c

H?

p

r

c cX 3200?

log(pn) + log 4
m

(31)

r

(32)

max{kr⇥ f (⇤? , ⇥? )k1 , kr⇤ f (⇤? , ⇥? )k1 } 

(33)

p

log(pn) + log 4.
mThen by Lemma 1 and the minimum bound on m we have that= c cX 3200?

with probability greater than 1 c1 exp( c2 m 2 ); we proceed with the proof conditioned on this event. Next,
note by our choice of we have that  ↵ /8 and thus the first half of the third condition for Lemma 2 holds.
It remains to show that R( )  ↵ /8. By our minimum bound on m and our choice of we have that
✓
◆
⇢
8
1
1
1
1+

min
,
(34)
↵
2cH ?
3c⌃? d 412c4⌃? c2⇥? c2X d2
and thus Lemma 4 applies, which gives
k k1  2cH ?

✓

1+

8
a

◆

.

(35)

Appendices to “Sparse Gaussian conditional random fields”

Therefore, the assumption of k k1  d1 min{ 3c1⌃? , c⇥2? } holds and we apply Lemma 3 to establishkR( )k1  206c4⌃? c2⇥? c2X d2 k k21
✓
◆2
8
4
2
2 2
2
?
 824c⌃? c⇥? cX d cH 1 +
a
"
✓
◆2 #
8
↵
4
2
2 2
 824c⌃? c⇥? cX d cH ? 1 +
a
8


(36)

↵
.
8

˜ ⇥)
˜ is unique and thus combined with the above
Finally, note that our Assumption 2 implies that the solution (⇤,˜ ⇥)˜ = (⇤,
ˆ ⇥)ˆ andderivation and Assumption 3, we the conditions for Lemma 2 and thus we conclude that (⇤,the thus claim 1 and 2 are satisfied.
C.1. Proofs of Lemmas
Proof. (of Lemma 1) Let X be given and assuming that Y is generated according to our model. We first consider
P (kr⇥ f (⇤? , ⇥? )k1 > ✏); as shown in Appendix A, we haver⇥ f (⇤, ⇥) = 2Sxy + 2Sxx ⇥⇤ 1 .
Writing
thus

⇤

=

(37)

⇥⇤ ⇤⇤ 1 and ⌃⇤ = ⇤⇤ 1 , we have Y = X ⇤ + Z where Z 2 Rm⇥p has rows Zi ⇠ N (0, ⌃⇤ ), and222Sxy + 2Sxx ⇥⇤ ⇤⇤ 1 = (X T Y X T X ⇤ ) = X T Z.
(38)
m
m
p
By our assumptions that kXj k2 / n < cX for all columns of X and the maximum diagonal entry ⌃? is c2 ? , we
have
✓
◆
2 T
4c2 ? c2X
Var
Xi Z j (39)mm
for any columns Xi and Zj . Therefore by the union bound and Gaussian tail probability we have
✓
◆
⇢
1
m✏2
P k X T Zk1 > ✏  2np expm
8c2 ⇤ c2X

(40)

Next, we consider P (kr⇤ f (⇤? , ⇥? )k1 > ✏) and again from Appendix A, we haver⇤ f (⇤, ⇥) = Syy⇤ 1⇤ 1 ⇥T Sxx ⇥⇤ 1(41)

which we can rewrite as

1 T
Z Z ⌃⇤ .
m
Now we can apply Lemma 1 in (Ravikumar et al., 2011) and arrive at the desired bound⇢1m✏2
P (k Z T Z ⌃⇤ k1 > ✏) < 4p2 expm3200c2 ?
Syy⇤⇤ 1⇤⇤ 1 ⇥⇤T Sxx ⇥⇤ ⇤⇤ 1 =

(42)

(43)

for 0 < ✏ < 40c ? .
Proof. (of Lemma 2)
The proof here proceeds exactly as in (Wainwright, 2009) and (Ravikumar et al., 2011), so we describe it relatively
quickly. The goal is to show that when solve the restricted problem, the resulting ✓˜ (which is zero outside thesupport set S) is also optimal for the full `1 problem. Defining = ✓˜ ✓? , the the full `1 optimization problem
can be written as
r2✓ f (✓? ) + r✓ f (✓? ) R( ) + z = 0.
(44)

Appendices to “Sparse Gaussian conditional random fields”

ˆ Furthermore,
If we can show that kzk1 < 1, then ✓˜ is an optimal solution to the original `1 problem, so ✓˜ = ✓.
?
ˆthe solution ✓ cannot have support outside the support of ✓ .
We can write the above optimality condition in terms of S and S̄, using H = r2✓ f (✓? ) and g = r✓ f (✓? ) for
simplicity





HSS HS S̄
gS
R( )S
zS
S
+
+
+
(45)
HS̄S HS̄ S̄
0gS̄
R( )S̄zS̄
Using the fact that1
S = HSS (R(

)S

gS

zS )

(46)

we can solve forzS̄ gives
zS̄ =

1

HS̄S

S +

1

(R( )S

Thus
kzS̄ k1 

2gS̄ ) =↵

1

HS̄S HSS1 (gS

R( )S ) + HS̄S HSS1 zS +(kgk1 + kR( )k1 ) + 1

↵

2

↵↵
+1
4

1

(R( )SgS )↵ < 1.

(47)(48)

Proof. (of Lemma 3) Since R( ) is the residual of the first order Taylor expansion of the likelihood gradient, by
the mean value theorem we have that the exists t 2 [0, 1] such that
R⇤ (

⇤,

2
?
⇥ ) = d (r⇤ f (⇤ + t

⇤, ⇥

?

+t

⇤ );

⇤,

⇥)

(49)

and similarly for R⇥ ( ⇤ , ⇥ ). The first and second di↵erentials of these gradient terms are given by (note thatsince these are the di↵erentials of a matrix-valued function, we cannot simplify as many of the expressions as we
did for the di↵erential of the likelihood function)d(r⇤ f (⇤? , ⇥? ); U, V ) = ⇤ 1 U ⇤ 1 + ⇤ 1 U ⇤ 1 ⇥T Sxx ⇥⇤ 1 + ⇤ 1 ⇥T Sxx ⇥⇤ 1 U ⇤ 1
d2 (r⇤ f (⇤? , ⇥? ); U, V ) =⇤ 1 V T Sxx ⇥⇤ 1

⇤ 1 ⇥T Sxx V ⇤ 1

2⇤ 1 U ⇤ 1 U ⇤ 1

2⇤ 1 U ⇤ 1 U ⇤ 1 ⇥T Sxx ⇥⇤ 1

2⇤ 1 U ⇤ 1 ⇥T Sxx ⇥⇤ 1 U ⇤ 1

2⇤ 1 ⇥T Sxx ⇥⇤ 1 U ⇤ 1 U ⇤ 1 +2⇤ 1 U ⇤ 1 V T Sxx ⇥⇤ 1 + 2⇤ 1 U ⇤ 1 ⇥T Sxx V ⇤ 1 +?

?

d(r⇥ f (⇤ , ⇥ ); U, V ) =

2⇤

1

2⇤

1

T

1

T

1

V Sxx ⇥⇤V Sxx V ⇤

U⇤

1

+ 2⇤

1

T

⇥ Sxx V ⇤

1

U⇤

(50)

1

2Sxx ⇥⇤ 1 U ⇤ 1 + 2Sxx V ⇤ 1d2 (r⇥ f (⇤? , ⇥? ); U, V ) = 4Sxx ⇥⇤ 1 U ⇤ 1 U ⇤ 1

4Sxx V ⇤ 1 U ⇤ 1

For example, R⇤ ( ⇤ , ⇥ ) is equal to the second expression with the ⇤? terms replaced by ⇤? + ⇤ , the ⇥?
terms replaced by ⇥? + ⇥ and U and V replaced by ⇤ and ⇥ respectively. To bound R( ), we bound eachof these terms individually.
Although the expression is rather lengthy, note that each term in the second di↵erentials has a quadratic expression in ⇤ and ✓ , with at most four (⇤? + t ⇤ ) 1 terms, two ⇥? + t ⇥ terms and one Sxx term. Furthermore,
we use the fact that
kABCk1 = k(C T ⌦ A) vec(B)k1  |||C|||1 |||A|||1 kBk1
(51)
to place the vector infinity norm around the Sxx term in all cases, since kSxx k1  c2X . Thus, each term in thesecond di↵erential is bounded by
c2X |||(⇤? + t
Now, first note that since⇤)1 4|||1 |||⇥? + t

2⇥ |||1 |||

|||1

(52)

as a most d entries per column||| |||1  dk k1 .

(53)

Appendices to “Sparse Gaussian conditional random fields”

Now, note that(⇤? + t

⇤)

1

⇤)

1

= (I + t⇤? 1

and(I + t⇤

? 1

=

1
X

⇤)

1

( 1)i (t⇤? 1

(54)

⇤)

i

(55)

i=1

so that|||(⇤? + t⇤)

1|||1  |||⇤? 1 |||1


1

1
X
i=1

i

|||⇤? 1 |||1 |||

i

c ⌃?
3c⌃?

.
c⌃? dk k1
2

Furthermore,|||⇥? + t

⇤ |||1 k?
⇥ |||1  |||⇥ |||1 + |||⇥ |||1  c⇥? +1
3c⇥?
dk k1 
.
2
2

(56)

(57)

Combining these expressions results in the boundkR( )k1  206c4⌃? c2⇥? c2X d2 k k21(58)

as required.
˜ ⇥)˜ be the solution to the
Proof. (of Lemma 4) Let (⇤? , ⇥? ) be the true parameters with support S and (⇤,optimization problem restricted to this support set. Our goal is to bound k k1 where= [ ⇤ ⇥ ] with
?
?
˜
˜
=
⇤
⇤
and
=
⇥
⇥
.
⇤
⇥
Definer := 2cH ? (max{kr⇤ f (⇤? , ⇥? )k1 , kr⇥ f (⇤? , ⇥? )k1 } + )(59)

and note that by assumption we haver  2cH ? (minTo bound k k1 observe that we have⇢13c⌃? d,1
412c4⌃? c2⇥? cX d2

)

(60)

C = 0 and? 1
S = HSS (RS () + GSZS )(61)

as shown in the proof for Lemma 2. Our approach will be the same as that of (Ravikumar et al., 2011), using
Brouwer’s fixed point theorem. To do so, note that we can view the RHS of the above equation as a continuousfunction of and thus by Brouwer’s fixed point theorem on a compact set, it suffices to show that if show that⇤ 1
if k S k1  r then kHSS(RS GSZS )k1  r as this implies that there is a solution to this equation such˜ ⇥).
˜that k S k  r and by uniqueness (from Assumption 2) this solution must be (⇤,Taking infinity norm, we havek⇤ 1
S k1  kHSS k1 kR(⇤ 1)k1 + kHSS
k1 kGS

ZS k1

(62)

For the first term, through application of the bound on R( ) and by assumption on k k1
⇤ 1
kHSS
k1 kR( )k1  206H ⇤ c4⌃? c2⇥? cX d2 k k21 

And for the second term
and thus the claim is proven.

⇤ 1
kHSS
k1 kGSZS k1  H ⇤ (kGk1 + ) 

r
2

r
2

(63)

(64)

Appendices to “Sparse Gaussian conditional random fields”

y1

y2

y3

x1

x2

x3

Figure 1. The chain CRF with 3 input variables and 3 output variables.

1
0.9

Mutual incoherence

0.8
0.7
0.6
0.5
0.4
0.3

β = 0.1β=1
β=5β = 10

0.2
0.1
0

0

0.05

0.1

0.15
ρ

0.2

0.25

0.3

Figure 2. The mutual incoherence condition |||HS̄S (HSS ) 1 |||1 while varying ⇢ and .

D. Mutual Incoherence for the Chain CRF
In this section we consider the mutual incoherence condition for the chain CRF, illustrated in Figure 1. For
simplicity, we consider a class of models parameterized by two variables: ⇢ describing the conditional dependence
between the output variables and describing the relative influence of the input variables on the output variables.
In particular, the class of models specified by
2
3
2
3
1 ⇢ 0
⇢
0
0
05
⇤? = 4⇢ 1 ⇢5 ⇥? = 4 0 ⇢
(65)
0 ⇢ 1
0
0 ⇢
with positive ⇢ and .

We are interested in characterizing the range over which the mutual incoherence condition|||HS̄S (HSS ) 1 |||1 < 1(66)

holds. Note that the Hessian (given in Appendix A) depends not only on these parameters, but also on theempirical covariance of the input features, Sxx ; for the purpose of this illustration, we take Sxx to be the identity
matrix, representing an ideal case in which the input features are perfectly uncorrelated. Under these conditions,
we can see from Figure 2 that mutual incoherence indeed holds over a range of the parameters. However, as ⇢increases and the output variables become more correlated, we approach the boundary at which this assumption
is no longer valid.

Appendices to “Sparse Gaussian conditional random fields”

References
Hsieh, C.-J., Sustik, M. A., Dhillon, I. S., and Ravikumar, P. Sparse inverse covariace matrix estimation using
quadratic approximation. In Neural Information Processing Systems, 2011.
Magnus, X and Neudecker, Heinz. Matrix di↵erential calculus. New York, 1988.
Ravikumar, Pradeep, Wainwright, Martin J, Raskutti, Garvesh, and Yu, Bin. High-dimensional covarianceestimation by minimizing 1-penalized log-determinant divergence. Electronic Journal of Statistics, 5:935–980,
2011.
Wainwright, M.J. Sharp thresholds for high-dimensional and noisy sparsity recovery using `1 -constrained
quadratic programming (Lasso). IEEE Transactions on Information Theory, 55(5):2183–2202, 2009.