The Multiple Quantile Graphical Model
Alnur Ali
Machine Learning Department
Carnegie Mellon University
alnurali@cmu.edu

J. Zico Kolter
Computer Science Department
Carnegie Mellon University
zkolter@cs.cmu.edu

Ryan J. Tibshirani
Department of Statistics
Carnegie Mellon University
ryantibs@cmu.edu

Abstract
We introduce the Multiple Quantile Graphical Model (MQGM), which extends
the neighborhood selection approach of Meinshausen and Bühlmann for learning
sparse graphical models. The latter is defined by the basic subproblem of modeling the conditional mean of one variable as a sparse function of all others. Our
approach models a set of conditional quantiles of one variable as a sparse functionof all others, and hence offers a much richer, more expressive class of conditional
distribution estimates. We establish that, under suitable regularity conditions, the
MQGM identifies the exact conditional independencies with probability tending to
one as the problem size grows, even outside of the usual homoskedastic Gaussian
data model. We develop an efficient algorithm for fitting the MQGM using thealternating direction method of multipliers. We also describe a strategy for sampling from the joint distribution that underlies the MQGM estimate. Lastly, we
present detailed experiments that demonstrate the flexibility and effectiveness of
the MQGM in modeling hetereoskedastic non-Gaussian data.

1

Introduction

We consider modeling the joint distribution Pr(y1 , . . . , yd ) of d random variables, given n independent draws from this distribution y (1) , . . . , y (n) 2 Rd , where possibly dn. Later, we generalize
this setup and consider modeling the conditional distribution Pr(y1 , . . . , yd |x1 , . . . , xp ), given n
independent pairs (x(1) , y (1) ), . . . , (x(n) , y (n) ) 2 Rp+d . Our starting point is the neighborhood selection method [28], which is typically considered in the context of multivariate Gaussian data, and seen
as a tool for covariance selection [8]: when Pr(y1 , . . . , yd ) is a multivariate Gaussian distribution, it
is a well-known fact yj and yk are conditionally independent given the remaining variables if and ony
if the coefficent corresponding to yk is zero in the (linear) regression of yj on all other variables (e.g.,[22]). Therefore, in neighborhood selection we compute, for each k = 1, . . . , d, a lasso regression —in order to obtain a small set of conditional dependencies — of yk on the remaining variables, i.e.,
◆2
n ✓
X
X
(i)
(i)
minimizeyk✓kj yj+ k✓k k1 ,(1)
✓k 2Rdfor a tuning parameter

i=1

j6=k

> 0. This strategy can be seen as a pseudolikelihood approximation [4],Pr(y1 , . . . , yd ) ⇡

d
Y

k=1

Pr(yk |y¬k ),(2)where y¬k denotes all variables except yk . Under the multivariate Gaussian model for Pr(y1 , . . . , yd ),the conditional distributions Pr(yk |y¬k ), k = 1, . . . , d here are (univariate) Gaussians, and maximizing the pseudolikelihood in (2) is equivalent to separately maximizing the conditionals, as is precisely
done in (1) (with induced sparsity), for k = 1, . . . , d.
29th Conference on Neural Information Processing Systems (NIPS 2016), Barcelona, Spain.

Following the pseudolikelihood-based approach traditionally means carrying out three steps: (i) we
write down a suitable family of joint distributions for Pr(y1 , . . . , yd ), (ii) we derive the conditionals
Pr(yk |y¬k ), k = 1, . . . , d, and then (iii) we maximize each conditional likelihood by (freely) fitting
the parameters. Neighborhood selection, and a number of related approaches that came after it (seeSection 2.1), can be all thought of in this workflow. In many ways, step (ii) acts as the bottleneck
here, and to derive the conditionals, we are usually limited to a homoskedastic and parameteric family
for the joint distribution.
The approach we take in this paper differs somewhat substantially, as we begin by directly modelingthe conditionals in (2), without any preconceived model for the joint distribution — in this sense, it
may be seen a type of dependency network [13] for continuous data. We also employ heteroskedastic,
nonparametric models for the conditional distributions, which allows us great flexibility in learningthese conditional relationships. Our method, called the Multiple Quantile Graphical Model (MQGM),
is a marriage of ideas in high-dimensional, nonparametric, multiple quantile regression with those in
the dependency network literature (the latter is typically focused on discrete, not continuous, data).
An outline for this paper is as follows. Section 2 reviews background material, and Section 3 develops
the MQGM estimator. Section 4 studies basic properties of the MQGM, and establishes a structure
recovery result under appropriate regularity conditions, even for heteroskedastic, non-Gaussian data.
Section 5 describes an efficient ADMM algorithm for estimation, and Section 6 presents empiricalexamples comparing the MQGM versus common alternatives. Section 7 concludes with a discussion.

2

Background

2.1

Neighborhood selection and related methods

Neighborhood selection has motivated a number of methods for learning sparse graphical models. The
literature here is vast; we do not claim to give a complete treatment, but just mention some relevant
approaches. Many pseudolikelihood approaches have been proposed, see e.g., [35, 33, 12, 24, 17, 1].
These works exploit the connection between estimating a sparse inverse covariance matrix and regression, and they vary in terms of the optimization algorithms they use and the theoretical guarantees
they offer. In a clearly related but distinct line of research, [45, 2, 11, 36] proposed `1 -penalized
likelihood estimation in the Gaussian graphical model, a method now generally termed the graphicallasso (GLasso). Following this, several recent papers have extended the GLasso in various ways. [10]examined a modification based on the multivariate Student t-distribution, for robust graphical modeling. [37, 46, 42] considered conditional distributions of the form Pr(y1 , . . . , yd |x1 , . . . , xp ). [23]proposed a model for mixed (both continuous and discrete) data types, generalizing both GLasso andpairwise Markov random fields. [25, 26] used copulas for learning non-Gaussian graphical models.
A strength of neighborhood-based (i.e., pseudolikelihood-based) approaches lies in their simplicity;
because they essentially reduce to a collection of univariate probability models, they are in a sensemuch easier to study outside of the typical homoskedastic, Gaussian data setting. [14, 43, 44] elegantly studied the implications of using univariate exponential family models for the conditionals in(2). Closely related to pseudoliklihood approaches are dependency networks [13]. Both frameworks
focus on the conditional distributions of one variable given all the rest; the difference lies in whether
or not the model for conditionals stems from first specifying some family of joint distributions (pseudolikelihood methods), or not (dependency networks). Dependency networks have been thoroughly
studied for discrete data, e.g., [13, 29]. For continuous data, [40] proposed modeling the mean in a
Gaussian neighborhood regression as a nonparametric, additive function of the remaining variables,
yielding flexible relationships — this is a type of dependency network for continuous data (though it
is not described by the authors in this way). Our method, the MQGM, also deals with continuous
data, and is the first to our knowledge that allows for fully nonparametric conditional distributions, aswell as nonparametric contributions of the neighborhood variables, in each local model.
2.2

Quantile regression

In linear regression, we estimate the conditional mean of y|x1 , . . . , xp from samples. Similarly, in ↵quantile regression [20], we estimate the conditional ↵-quantile of y|x1 , . . . , xp for a given ↵ 2 [0, 1],formally Qy|x1 ,...,xp (↵). . , xp ) ↵}, by solving the convex optimization
Pn= inf{t :(i)Pr(yPp t|x1 , .(i)problem: minimize✓1)z} is the ↵i=1 ↵ (yj=1 ✓j xj ), where ↵ (z) = max{↵z, (↵2quantile loss (also called the “pinball” or “tilted absolute” loss). Quantile regression can be usefulwhen the conditional distribution in question is suspected to be heteroskedastic and/or non-Gaussian,e.g., heavy-tailed, or if we wish to understand properties of the distribution other than the mean,e.g., tail behavior. In multiple quantile regression, we solve several quantile regression problems
simultaneously, each corresponding to a different quantile level; these problems can be coupledsomehow to increase efficiency in estimation (see details in the next section). Again, the literature
on quantile regression is quite vast (especially that from econometrics), and we only give a shortreview here. A standard text is [18]. Nonparametric modeling of quantiles is a natural extension fromthe (linear) quantile regression approach outlined above; in the univariate case (one conditioning
variable), [21] suggested a method using smoothing splines, and [38] described an approach using
kernels. More recently, [19] studied the multivariate nonparametric case (more than one conditioning
variable), using additive models. In the high-dimensional setting, where p is large, [3, 16, 9] studied`1 -penalized quantile regression and derived estimation and recovery theory for non-(sub-)Gaussiandata. We extend results in [9] to prove structure recovery guarantees for the MQGM (in Section 4.3).

3

The multiple quantile graphical model

Many choices can be made with regards to the final form of the MQGM, and to help in understanding
these options, we break down our presentation in parts. First fix some ordered set A = {↵1 , . . . , ↵r }of quantile levels, e.g., A = {0.05, 0.10, . . . , 0.95}. For each variable yk , and each level ↵` , wemodel the conditional ↵` -quantile given the other variables, using an additive expansion of the form:
Qyk |y¬k (↵` ) = b⇤`k +

d
X

⇤
f`kj
(yj ),

(3)

j6=k

⇤
where b⇤`k 2 R is an intercept term, and f`kj
, j = 1, . . . , d are smooth, but not parametric in form. In
its most general form, the MQGM estimator is defined as a collection of optimization problems, overk = 1, . . . , d and ` = 1, . . . , r:
✓
◆ X⇣
n
⌘!
X
X
(i)
(i)
minimizeb`kf`kj (yj ) +. (4)↵ ` yk1 P1 (f`kj ) + 2 P2 (f`kj )b`k , f`kj 2F`kj ,i=1
j=1,...,d

j6=k

j6=k

Here 1 , 2 0 are tuning parameters, F`kj , j = 1, . . . , d are univariate function spaces, ! > 0 is
a fixed exponent, and P1 , P2 are sparsity and smoothness penalty functions, respectively. We give
three examples below; many other variants are also possible.

Example 1: basis expansion model Consider taking F`kj = span{ j1 , . . . , jm }, the span of mbasis functions, e.g., radial basis functions (RBFs) with centers placed at appropriate locations acrossthe domain of variable j, for each j = 1, . . . , d. This means that each f`kj 2 F`kj can be expressedTj
as f`kj (x) = ✓`kj(x), for a coefficient vector ✓`kj 2 Rm , where j (x) = ( j1 (x), . . . , jm (x)).
Also consider an exponent ! = 1, and the sparsity and smoothness penalties
P1 (f`kj ) = k✓`kj k2and

P2 (f`kj ) = k✓`kj k22 ,respectively, which are group lasso and ridge penalties, respectively. With these choices in place, the
MQGM problem in (4) can be rewritten in finite-dimensional form:
⇣
⌘ X⇣
⌘
2
minimize
b`k 1
✓`k +
(5)
↵` Yk
1 k✓`kj k2 + 2 k✓`kj k2 .
b`k , ✓`k =(✓`k1 ,...,✓`kd )

j6=k

Pn
Above, we have used the abbreviation ↵` (z) = i=1 ↵` (zi ) for a vector z = (z1 , . . . , zn ) 2 Rn ,(1)(n)
and also Yk = (yk , . . . , yk ) 2 Rn for the observations along variable k, 1 = (1, . . . , 1) 2 Rn , and
(i)
n⇥dm
2R
for the basis matrix, with blocks of columns to be understood as ij = (yj )T 2 Rm .
The basis expansion model is simple and tends to work well in practice, so we focus on it for most of
the paper. In principle, essentially all our results apply to the next two models we describe, as well.

Example 2: smoothing splines model Now consider taking F`kj = span{g1j , . . . , gnj }, the span(1)(n)
of m = n natural cubic splines with knots at yj , . . . , yj , for j = 1, . . . , d. As before, we canTjthen write f`kj (x) = ✓`kj g (x) with coefficients ✓`kj 2 Rn , for f`kj 2 F`kj . The work of [27], on
high-dimensional additive smoothing splines, suggests a choice of exponent ! = 1/2, and penalties
P1 (f`kj ) = kGj ✓`kj k22

T
and P2 (f`kj ) = ✓`kj⌦j ✓`kj ,3

for sparsity and smoothness, respectively, where Gj 2 Rn⇥n is a spline basis matrix with entries(i)
Gjii0 = gij0 (yj ), and ⌦j is the smoothing spline penalty matrix containing integrated products ofpairs of twice differentiated basis functions. The MQGM problem in (4) can be translated into a
finite-dimensional form, very similar to what we have done in (5), but we omit this for brevity.
Example 3: RKHS model Consider taking F`kj = Hj , a univariate reproducing kernel Hilbert
space (RKHS), with kernel function j (·, ·). The representer theorem allows
Pn us to express each(i)function f`kj 2 Hj in terms of the representers of evaluation, i.e., f`kj (x) = i=1 (✓`kj )i j (x, yj ),nfor a coefficient vector ✓`kj 2 R . The work of [34], on high-dimensional additive RKHS modeling,suggests a choice of exponent ! = 1, and sparsity and smoothness penaltiesq
T Kj✓
P1 (f`kj ) = kK j ✓`kj k2 and P2 (f`kj ) = ✓`kj
`kj ,
(i)

(i0 )

respectively, where K j 2 Rn⇥n is the kernel matrix with entries Kiij 0 = j (yj , yj ). Again, the
MQGM problem in (4) can be written in finite-dimensional form, now an SDP, omitted for brevity.

Structural constraints
Several structural constraints can be placed on top of the MQGM optimization problem in order to guide the estimated component functions to meet particular shaperequirements. An important example are non-crossing constraints (commonplace in nonparametric,
multiple quantile regression [18, 38]): here, we optimize (4) jointly over ` = 1, . . . , r, subject to
X
X
(i)
(i)
b`k +
f`kj (yj )  b`0 k +f`0 kj (yj ), for all ↵` < ↵`0 , and i = 1, . . . , n.
(6)
j6=k

j6=k

This ensures that the estimated quantiles obey the proper ordering, at the observations. For concreteness, we consider the implications for the basis regression model, in Example 1 (similar statements
hold for the other two models). For each ` = 1, . . . , r, denote by F`k (b`k , ✓`k ) the criterion in (5).
Introducing the non-crossing constraints requires coupling (5) over ` = 1, . . . , r, so that we now havethe following optimization problems, for each target variable k = 1, . . . , d:rXminimizeF`k (b`k , ✓`k ) subject to (1BkT + ⇥k )DT 0,(7)
Bk ,⇥k`=1

where we denote Bk = (b1k , . . . , brk ) 2 Rr , 2 Rn⇥dm the basis matrix as before, ⇥k 2 Rdm⇥rgiven by column-stacking ✓`k 2 Rdm , ` = 1, . . . , r, and D 2 R(r 1)⇥r is the usual discrete
difference operator. (The inequality in (7) is to be interpreted componentwise.) Computationally,
coupling the subproblems across ` = 1, . . . , r clearly adds to the overall difficulty of the MQGM, but
statistically this coupling acts as a regularizer, by constraining the parameter space in a useful way,thus increasing our efficiency in fitting multiple quantile levels from the given data.
(i)

(i0 )

For a triplet`, k, j, monotonicity constraints are also easy to add, i.e., f`kj (yj )  f`kj (yj ) for all
(i)
(i0 )
yj < yj . Convexity constraints, where we require f`kj to be convex over the observations, for aparticular `, k, j, are also straightforward. Lastly, strong non-crossing constraints, where we enforce(6) over all z 2 Rd (not just over the observations) are also possible with positive basis functions.

Exogenous variables and conditional random fields So far, we have considered modeling thejoint distribution Pr(y1 , . . . , yd ), corresponding to learning a Markov random field (MRF). It is not
hard to extend our framework to model the conditional distribution Pr(y1 , . . . , yd |x1 , . . . , xp ) given
some exogenous variables x1 , . . . , xp , corresponding to learning a conditional random field (CRF).
x
To extend the basis regression model, we introduce the additional parameters ✓`k
2 Rp in (5), and then⇥qTxloss now becomes ↵` (Yk b`k 1✓`k X✓`k ), where X 2 R
is filled with the exogenous
observations x(1) , . . . , x(n) 2 Rq ; the other models are changed similarly.

4

Basic properties and theory

4.1

Quantiles and conditional independence⇤
In the model (3), when a particular variable yj has no contribution, i.e., satisfied f`kj
= 0 across allquantile levels ↵` , ` = 1, . . . , r, what does this imply about the conditional independence between yk
and yj , given the rest? Outside of the multivariate normal model (where the feature transformations
need only be linear), nothing can be said in generality. But we argue that conditional independence canbe understood in a certain approximate sense (i.e., in a projected approximation of the data generating
model). We begin with a simple lemma. Its proof is elementary, and given in the supplement.

4

Lemma 4.1. Let U, V, W be random variables, and suppose that all conditional quantiles of U |V, W
do not depend on V , i.e., QU |V,W (↵) = QU |W (↵) for all ↵ 2 [0, 1]. Then U and V are conditionallyindependent given W .
By the lemma, if we knew that QU |V,W (↵) = h(↵, U, W ) for a function h, then it would follow that
U, V are conditionally independent given W (n.b., the converse is true, as well). The MQGM problemin (4), with sparsity imposed on the coefficients, essentially aims to achieve such a representation
for the conditional quantiles; of course we cannot use a fully nonparametric representation of theconditional distribution yk |y¬k and instead we use an r-step approximation to the conditional cumulative distribution function (CDF) of yk |y¬k (corresponding to estimating r conditional quantiles),and (say) in the basis regression model, limit the dependence on conditioning variables to be in termsof an additive function of RBFs in yj , j 6= k. Thus, if at the solution in (5) we find that ✓ˆkj` = 0,
` = 1, . . . , r, we may interpret this to mean that yk and yj are conditionally independent given theremaining variables, but according to the distribution defined by the projection of yk |y¬k onto thespace of models considered in (5) (r-step conditional CDFs, which are additive expansions in yj ,j 6= k). This interpretation is no more tenuous (arguably, less so, as the model space here is muchlarger) than that needed when applying standard neighborhood selection to non-Gaussian data.
4.2

Gibbs sampling and the “joint” distribution

When specifying a form for the conditional distributions in a pseudolikelihood approximation as in(2), it is natural to ask: what is the corresponding joint distribution? Unfortunately, for a generalcollection of conditional distributions, there need not exist a compatible joint distribution, evenwhen all conditionals are continuous [41]. Still, pseudolikelihood approximations (a special caseof composite likelihood approximations), possess solid theoretical backing, in that maximizing thepseudolikelihood relates closely to minimizing a certain (expected composite) Kullback-Leibler
divergence, measured to the true conditionals [39]. Recently, [7, 44] made nice progress in describing
specific conditions on conditional distributions that give rise to a valid joint distribution, though theirwork was specific to exponential families. A practical answer to the question of this subsection is to
use Gibbs sampling, which attempts to draw samples consistent with the fitted conditionals; this isprecisely the observation of [13], who show that Gibbs sampling from discrete conditionals converges
to a unique stationary distribution, although this distribution may not actually be compatible with theconditionals. The following result establishes the analogous claim for continuous conditionals; its
proof is in the supplement. We demonstrate the practical value of Gibbs sampling through various
examples in Section 6.
Lemma 4.2. Assume that the conditional distributions Pr(yk |y¬k ), k = 1, . . . , d take only positive
values on their domain. Then, for any given ordering of the variables, Gibbs sampling converges to aunique stationary distribution that can be reached from any initial point. (This stationary distribution
depends on the ordering.)
4.3

Graph structure recovery

When log d = O(n2/21 ), and we assume somewhat standard regularity conditions (listed as A1–A4
in the supplement), the MQGM estimate recovers the underlying conditional independencies withhigh probability (interpreted in the projected model space, as explained in Section 4.1). Importantly,
we do not require a Gaussian, sub-Gaussian, or even parametric assumption on the data generatingprocess; instead, we assume i.i.d. draws y (1) , . . . , y (n) 2 Rd , where the conditional distributions
yk |y¬k have quantiles specified by the model in (3) for k = 1, . . . , d, ` = 1, . . . , r, and further, each
⇤
T
j
⇤
f`kj
(x) = ✓`kj(x)⇤ for coefficients ✓`kj
2 Rm , j = 1, . . . , d, as in the basis expansion model.
Let E ⇤ denote the corresponding edge set of conditional dependencies from these neighborhood⇤⇤models, i.e., {k, j} 2 E ⇤ () max`=1,...,r max{k✓`kj
k2 , |✓`jkk2 } > 0. We define the estimatededge set Ê in the analogous way, based on the solutionin(5).
Without
a loss of generality, we assumep
the features have been scaled to satisfy k j k  n for all j = 1, . . . , dm. The following is ourrecovery result; its proof is provided in the supplement.

Theorem 4.3. Assume log d = O(n2/21 ), and conditions A1–A4 in the supplement. Assume that⇤the tuning parameters 1 , 2 satisfy 1 ⇣ (mn log(d2 mr/ ) log3 n)1/2 and 2 = o(n41/42 /✓max),⇤⇤
where ✓max = max`,k,j k✓`kj k2 . Then for n sufficiently large, the MQGM estimate in (5) exactly
recovers the underlying conditional dependencies, i.e., Ê = E ⇤ , with probability at least 1.
5

The theorem shows that the nonzero pattern in the MQGM estimate identifies, with high probability,
the underlying conditional independencies. But to be clear, we emphasize that the MQGM estimate
is not an estimate of the inverse covariance matrix itself (this is also true of neighborhood regression,SpaceJam of [40], and many other methods for learning graphical models).

5

Computational approach

By design, the MQGM problem in (5) separates into d subproblems, across k = 1, . . . , d (it therefore
suffices to consider only a single subproblem, so we omit notational dependence on k for auxiliary variables). While these subproblems are challenging for off-the-shelf solvers (even for only
moderately-sized graphs), the key terms here all admit efficient proximal operators [32], which makes
operator splitting methods like the alternating direction method of multipliers [5] a natural choice.
As an illustration, we consider the non-crossing constraints in the basis regression model below.
Reparameterizing our problem, so that we may apply ADMM, yields:
Pr Pd
2
T
2
minimize
A (Z) + 1`=1
j=1 kW`j k2 + 2 kW kF + I+ (V D )⇥k ,Bk ,V,W,Z(8)
subject to
V = 1BkT + ⇥k , W = ⇥k , Z = Yk 1T 1BkT
⇥k ,
Pr Pd
where for brevity A (A) = `=1 j=1 ↵` (A`j ), and I+ (·) is the indicator function of the spaceof elementwise nonnegative matrices. The augmented Lagrangian associated with (8) is:L⇢ (⇥k , Bk , V, W, Z, UV , UW , UZ ) =⇢⇣+k1BkT + ⇥k
2

A (Z) +

1r X
dX`=1 j=1

kW`j k2 +2

2

kW k2F + I+ (V DT )⌘Z + UZ k2F ,(9)where ⇢ > 0 is the augmented Lagrangian parameter, and UV , UW , UZ are dual variables corresponding to the equality constraints on V, W, Z, respectively. Minimizing (9) over V yields:
V + UV k2F + k⇥kV

W + UW k2F + kYk 1T

1BkT

⇥k

Piso 1BkT + ⇥k + UV ,(10)where Piso (·) denotes the row-wise projection operator onto the isotonic cone (the space of componentwise nondecreasing vectors), an O(nr) operation here [15]. Minimizing (9) over W`j yields the
update:✓◆(⇥k )`j + (UW )`j
1 /⇢
W`j
1
,
(11)
1 + 2 /⇢k(⇥k )`j + (UW )`j k2 +where (·)+ is the positive part operator. This can be seen by deriving the proximal operator of thefunction f (x) = 1 kxk2 + ( 2 /2)kxk22 . Minimizing (9) over Z yields the update:Z

prox(1/⇢) A (Yk 1T

1bTk

(12)

⇥k + UZ ),where proxf (·) denotes the proximal operator of a function f . For the multiple quantile loss function
A , this is a kind of generalized soft-thresholding. The proof is given in the supplement.
Lemma 5.1. Let P+ (·) and P (·) be the elementwise positive and negative part operators, respectively, and let a = (↵1 , . . . , ↵r ). Then proxt A (A) = P+ (A t1aT ) + P (A t1aT ).
Finally, differentiation in (9) with respect to Bk and ⇥k yields the simultaneous updates:

 T
1✓
T
1
⇥k
+ 12 I1T[I 0] (W UW ) +
BkT
1T
1T 1
2
[

T

1] (Yk 1T

Z + UZ + V◆
UV ) .

(13)

A complete description of our ADMM algorithm for solving the MQGM problem is in the supplement.
Gibbs sampling Having fit the conditionals yk |y¬k , k = 1, . . . d, we may want to make predictions
or extract joint distributions over subsets of variables. As discussed in Section 4.2, there is no general
analytic form for these joint distributions, but the pseudolikelihood approximation underlying the
MQGM suggests a natural Gibbs sampler. A careful implementation that respects the additive modelin (3) yields a highly efficient Gibbs sampler, especially for CRFs; the supplement gives details.
6

6

Empirical examples

6.1

Synthetic data

We consider synthetic examples, comparing the MQGM to neighborhood selection (MB), the graphical lasso (GLasso), SpaceJam [40], the nonparanormal skeptic [26], TIGER [24], and neighborhood
selection using the absolute loss (Laplace).
Ring example As a simple but telling example, we drew n = 400 samples from a “ring” distribution
in d = 4 dimensions. Data were generated by drawing a random angle ⌫ ⇠ Uniform(0, 1), a random
radius R ⇠ N (0, 0.1), and then computing the coordinates y1 = R cos ⌫, y2 = R sin ⌫ andy3 , y4 ⇠ N (0, 1), i.e., y1 and y2 are the only dependent variables here. The MQGM was used with
m = 10 basis functions (RBFs), and r = 20 quantile levels. The left panel of Figure 1 plots samples(blue) of the coordinates y1 , y2 as well as new samples from the MQGM (red) fitted to these same
(blue) samples, obtained by using our Gibbs sampler; the samples from the MQGM appear to closely
match the samples from the underlying ring. The main panel of Figure 1 shows the conditional
dependencies recovered by the MQGM, SpaceJam, GLasso, and MB (plots for the other methods are
given in the supplement), when run on the ring data. We visualize these dependencies by forming ad ⇥ d matrix with the cell (j, k) set to black if j, k are conditionally dependent given the others, andwhite otherwise. Across a range of tuning parameters for each method, the MQGM is the only one
that successfully recovers the underlying conditional dependencies, at some point along its solution
path. In the supplement, we present an evaluation of the conditional CDFs given by each method,when run on the ring data; again, the MQGM performs best in this setting.
Larger examples To investigate performance at larger scales, we drew n 2 {50, 100, 300} samples
from a multivariate normal and Student t-distribution (with 3 degrees of freedom), both in d = 100
dimensions, both parameterized by a random, sparse, diagonally dominant d ⇥ d inverse covariancematrix, following the procedure in [33, 17, 31, 1]. Over the same set of sample sizes, with d = 100, we
also considered an autoregressive setup in which we drew samples of pairs of adjacent variables fromthe ring distribution. In all three data settings (normal, t, and autoregressive), we used m = 10 andr = 20 for the MQGM. To summarize the performances, we considered a range of tuning parameters
for each method, computed corresponding false and true positive rates (in detecting conditional
dependencies), and then computed the corresponding area under the curve (AUC), following, e.g.,[33, 17, 31, 1]. Table 1 reports the median AUCs (across 50 trials) for all three of these examples; theMQGM outperforms all other methods on the autoregressive example; on the normal and Student t
examples, it performs quite competitively.
Truth

λ1 =8.00000

λ1 =16.00000

λ1 =32.00000

λ1 =64.00000

λ1 =128.00000

7ruth

λ1 =0.00781

λ1 =0.01562

λ1 =0.03125

λ1 =0.06250

λ1 =0.12500

0.5
0.0

GLasso−1.5−1.5 −1.0 −0.5 0.0

0.5

1.0

λ1 =0.12500

λ1 =0.25000

λ1 =0.50000

λ1 =1.00000

λ1 =2.00000

Truth

λ1 =0.06250

λ1 =0.12500

λ1 =0.25000

λ1 =0.50000

λ1 =1.00000

6SaceJam

−0.5
−1.0

Truth

0B

04G0

1.5
1.0

1.5

Figure 1: Left: data from the ring distribution (blue) as well as new samples from the MQGM (red) fitted tothe same (blue) data, obtained by using our Gibbs sampler. Right: conditional dependencies recovered by the
MQGM, MB, GLasso, and SpaceJam on the ring data; black means conditional dependence. The MQGM is the
only method that successfully recovers the underlying conditional dependencies along its solution path.
Table 1: AUC values for the MQGM, MB, GLasso, SpaceJam, the nonparanormal skeptic, TIGER, andLaplace for the normal, t, and autoregressive data settings; higher is better, best in bold.
MQGM
MB
GLassoSpaceJam
Nonpara.
TIGER
Laplacen = 50
0.953
0.850
0.908
0.889
0.881
0.732
0.803

Normaln = 100
0.976
0.959
0.964
0.968
0.962
0.921
0.931

n = 300
0.988
0.994
0.998
0.997
0.996
0.996
0.989

Student t
n = 100
0.947
0.923
0.605
0.965
0.942
0.873
0.876

n = 50
0.928
0.844
0.691
0.893
0.862
0.420
0.800

7

n = 300
0.981
0.988
0.965
0.993
0.998
0.989
0.991

n = 50
0.726
0.532
0.541
0.624
0.545
0.503
0.530

Autoregressiven = 100n = 300
0.754
0.955
0.563
0.725
0.620
0.711
0.708
0.854
0.590
0.612
0.518
0.718
0.554
0.758

1

10

8

Figure 2: Top panel and bottom row,
middle panel: conditional dependencies recovered by the MQGM on the
flu data; each of the first ten cells corresponds to a region of the U.S., and blackmeans dependence. Bottom row, left
panel: wallclock time (in seconds) for
solving one subproblem using ADMM
versus SCS. Bottom row, right panel:samples from the fitted marginal distribution of the weekly flu incidence rates
at region 6; samples at larger quantiles
are shaded lighter, and the median is indarker blue.

2

5
7

3

9

4

6

1

10

8

2

5
7

3

9

4

6

1

6

11

16

1

2bMective vaOue

04G0
SCS
10

6

3

11

10

2

10

1

16

6.2

5egiRn 6

9

4

0

50

100
150
SecoQGs

200

250% Rf flu-like sympWRms

10

8
7
6
5
4
3
2
1
0
30

40

50

8

18

28

Week

Modeling flu epidemics

We study n = 937 weekly flu incidence reports from September 28, 1997 through August 30,
2015, across 10 regions in the United States (see the top panel of Figure 2), obtained from [6]. We
considered d = 20 variables: the first 10 encode the current week’s flu incidence (precisely, thepercentage of doctor’s visits in which flu-like symptoms are presented) in the 10 regions, and the last
10 encode the same but for the prior week. We set m = 5, r = 99, and also introduced exogenous
variables to encode the week numbers, so p = 1. Thus, learning the MQGM here corresponds
to learning the structure of a spatiotemporal graphical model, and reduces to solving 20 multiple
quantile regression subproblems, each of dimension (19 ⇥ 5 + 1) ⇥ 99 = 9504. All subproblems
took about 1 minute on a 6 core 3.3 Ghz Core i7 X980 processor.
The bottom left panel in Figure 2 plots the time (in seconds) taken for solving one subproblem using
ADMM versus SCS [30], a cone solver that has been advocated as a reasonable choice for a classof problems encapsulating (4); ADMM outperforms SCS by roughly two orders of magnitude. The
bottom middle panel of Figure 2 presents the conditional independencies recovered by the MQGM.
Nonzero entries in the upper left 10 ⇥ 10 submatrix correspond to dependencies between the yk
variables for k = 1, . . . , 10; e.g., the nonzero (0,2) entry suggests that region 1 and 3’s flu reports are
dependent. The lower right 10 ⇥ 10 submatrix corresponds to the yk variables for k = 11, . . . , 20,
and the nonzero banded entries suggest that at any region the previous week’s flu incidence (naturally)
influences the next week’s. The top panel of Figure 2 visualizes these relationships by drawing anedge between dependent regions; region 6 is highly connected, suggesting that it may be a bellwether
for other regions, roughly in keeping with the current understanding of flu dynamics. To draw samples
from the fitted distributions, we ran our Gibbs sampler over the year, generating 1000 total samples,making 5 passes over all coordinates between each sample, and with a burn-in period of 100 iterations.
The bottom right panel of Figure 2 plots samples from the marginal distribution of the percentagesof flu reports at region 6 (other regions are in the supplement) throughout the year, revealing theheteroskedastic nature of the data.
For space reasons, our last example, on wind power data, is presented in the supplement.

7

Discussion

We proposed and studied the Multiple Quantile Graphical Model (MQGM). We established theoreticaland empirical backing to the claim that the MQGM is capable of compactly representing relationships
between heteroskedastic non-Gaussian variables. We also developed efficient algorithms for both
estimation and sampling in the MQGM. All in all, we believe that our work represents a step forwardin the design of flexible yet tractable graphical models.
Acknowledgements AA was supported by DOE Computational Science Graduate Fellowship DEFG02-97ER25308. JZK was supported by an NSF Expeditions in Computation Award, CompSustNet,
CCF-1522054. RJT was supported by NSF Grants DMS-1309174 and DMS-1554123.
8

References[1] Alnur Ali, Kshitij Khare, Sang-Yun Oh, and Bala Rajaratnam. Generalized pseudolikelihood methods for inverse covariance estimation.
Technical report, 2016. Available at http://arxiv.org/pdf/1606.00033.pdf.
[2] Onureena Banerjee, Laurent El Ghaoui, and Alexandre d’Aspremont. Model selection through sparse maximum likelihood estimation
for multivariate Gaussian or binary data. Journal of Machine Learning Research, 9:485–516, 2008.
[3] Alexandre Belloni and Victor Chernozhukov. `1 -penalized quantile regression in high-dimensional sparse models. Annals of Statistics,
39(1):82–130, 2011.
[4] Julian Besag. Spatial interaction and the statistical analysis of lattice systems. Journal of the Royal Statistical Society: Series B, 36(2):
192–236, 1974.
[5] Stephen Boyd, Neal Parikh, Eric Chu, Borja Peleato, and Jonathan Eckstein. Distributed optimization and statistical learning via thealternating direction method of multipliers. Foundations and Trends in Machine Learning, 3(1):1–122, 2011.
[6] Centers for Disease Control and Prevention (CDC). Influenza national and regional level graphs and data, August 2015. URL http://gis.cdc.gov/grasp/fluview/fluportaldashboard.html.
[7] Shizhe Chen, Daniela Witten, and Ali Shojaie. Selection and estimation for mixed graphical models. Biometrika, 102(1):47–64, 2015.
[8] Arthur Dempster. Covariance selection. Biometrics, 28(1):157–175, 1972.
[9] Jianqing Fan, Yingying Fan, and Emre Barut. Adaptive robust variable selection. Annals of Statistics, 42(1):324–351, 2014.
[10] Michael Finegold and Mathias Drton. Robust graphical modeling of gene networks using classical and alternative t-distributions. Annalsof Applied Statistics, 5(2A):1057–1080, 2011.
[11] Jerome Friedman, Trevor Hastie, and Robert Tibshirani. Sparse inverse covariance estimation with the graphical lasso. Biostatistics, 9(3):432–441, 2008.
[12] Jerome Friedman, Trevor Hastie, and Robert Tibshirani. Applications of the lasso and grouped lasso to the estimation of sparse graphicalmodels. Technical report, 2010. Available at http://statweb.stanford.edu/~tibs/ftp/ggraph.pdf.
[13] David Heckerman, David Maxwell Chickering, David Meek, Robert Rounthwaite, and Carl Kadie. Dependency networks for inference,collaborative filtering, and data visualization. Journal of Machine Learning Research, 1:49–75, 2000.
[14] Holger Höfling and Robert Tibshirani. Estimation of sparse binary pairwise Markov networks using pseudo-likelihoods. Journal of
Machine Learning Research, 10:883–906, 2009.
[15] Nicholas Johnson. A dynamic programming algorithm for the fused lasso and `0 -segmentation. Journal of Computational and Graphical
Statistics, 22(2):246–260, 2013.
[16] Kengo Kato. Group lasso for high dimensional sparse quantile regression models. Technical report, 2011. Available at http://arxiv.
org/pdf/1103.1458.pdf.
[17] Kshitij Khare, Sang-Yun Oh, and Bala Rajaratnam. A convex pseudolikelihood framework for high dimensional partial correlation
estimation with convergence guarantees. Journal of the Royal Statistical Society: Series B, 77(4):803–825, 2014.
[18] Roger Koenker. Quantile Regression. Cambridge University Press, 2005.
[19] Roger Koenker. Additive models for quantile regression: Model selection and confidence bandaids. Brazilian Journal of Probability andStatistics, 25(3):239–262, 2011.
[20] Roger Koenker and Gilbert Bassett. Regression quantiles. Econometrica, 46(1):33–50, 1978.
[21] Roger Koenker, Pin Ng, and Stephen Portnoy. Quantile smoothing splines. Biometrika, 81(4):673–680, 1994.
[22] Steffen Lauritzen. Graphical models. Oxford University Press, 1996.
[23] Jason Lee and Trevor Hastie. Structure learning of mixed graphical models. In Proceedings of the 16th International Conference on
Artificial Intelligence and Statistics, pages 388–396, 2013.
[24] Han Liu and Lie Wang. TIGER: A tuning-insensitive approach for optimally estimating Gaussian graphical models. Technical report,
2012. Available at http://arxiv.org/pdf/1209.2437.pdf.
[25] Han Liu, John Lafferty, and Larry Wasserman. The nonparanormal: Semiparametric estimation of high dimensional undirected graphs.
Journal of Machine Learning Research, 10:2295–2328, 2009.
[26] Han Liu, Fang Han, Ming Yuan, John Lafferty, and Larry Wasserman. High-dimensional semiparametric Gaussian copula graphicalmodels. The Annals of Statistics, pages 2293–2326, 2012.
[27] Lukas Meier, Sara van de Geer, and Peter Buhlmann. High-dimensional additive modeling. Annals of Statistics, 37(6):3779–3821, 2009.
[28] Nicolai Meinshausen and Peter Bühlmann. High-dimensional graphs and variable selection with the lasso. Annals of Statistics, 34(3):1436–1462, 2006.
[29] Jennifer Neville and David Jensen. Dependency networks for relational data. In Proceedings of Fourth IEEE International Conference
on the Data Mining, pages 170–177. IEEE, 2004.
[30] Brendan O’Donoghue, Eric Chu, Neal Parikh, and Stephen Boyd. Operator splitting for conic optimization via homogeneous self-dual
embedding. Technical report, 2013. Available at https://stanford.edu/~boyd/papers/pdf/scs.pdf.
[31] Sang-Yun Oh, Onkar Dalal, Kshitij Khare, and Bala Rajaratnam. Optimization methods for sparse pseudolikelihood graphical model
selection. In Advances in Neural Information Processing Systems 27, pages 667–675, 2014.
[32] Neal Parikh and Stephen Boyd. Proximal algorithms. Foundations and Trends in Optimization, 1(3):123–231, 2013.
[33] Jie Peng, Pei Wang, Nengfeng Zhou, and Ji Zhu. Partial correlation estimation by joint sparse regression models. Journal of the American
Statistical Association, 104(486):735–746, 2009.
[34] Garvesh Raskutti, Martin Wainwright, and Bin Yu. Minimax-optimal rates for sparse additive models over kernel classes via convex
programming. Journal of Machine Learning Research, 13:389–427, 2012.
[35] Guilherme Rocha, Peng Zhao, and Bin Yu. A path following algorithm for sparse pseudo-likelihood inverse covariance estimation(SPLICE). Technical report, 2008. Available at https://www.stat.berkeley.edu/~binyu/ps/rocha.pseudo.pdf.
[36] Adam Rothman, Peter Bickel, Elizaveta Levina, and Ji Zhu. Sparse permutation invariant covariance estimation. Electronic Journal of
Statistics, 2:494–515, 2008.
[37] Kyung-Ah Sohn and Seyoung Kim. Joint estimation of structured sparsity and output structure in multiple-output regression via inverse
covariance regularization. In Proceedings of the 15th International Conference on Artificial Intelligence and Statistics, pages 1081–1089,
2012.
[38] Ichiro Takeuchi, Quoc Le, Timothy Sears, and Alexander Smola. Nonparametric quantile estimation. Journal of Machine Learning
Research, 7:1231–1264, 2006.
[39] Cristiano Varin and Paolo Vidoni. A note on composite likelihood inference and model selection. Biometrika, 92(3):519–528, 2005.
[40] Arend Voorman, Ali Shojaie, and Daniela Witten. Graph estimation with joint additive models. Biometrika, 101(1):85–101, 2014.
[41] Yuchung Wang and Edward Ip. Conditionally specified continuous distributions. Biometrika, 95(3):735–746, 2008.
[42] Matt Wytock and Zico Kolter. Sparse Gaussian conditional random fields: Algorithms, theory, and application to energy forecasting. In
Proceedings of the 30th International Conference on Machine Learning, pages 1265–1273, 2013.
[43] Eunho Yang, Pradeep Ravikumar, Genevera Allen, and Zhandong Liu. Graphical models via generalized linear models. In Advances in
Neural Information Processing Systems 25, pages 1358–1366, 2012.
[44] Eunho Yang, Pradeep Ravikumar, Genevera Allen, and Zhandong Liu. Graphical models via univariate exponential family distributions.
Journal of Machine Learning Research, 16:3813–3847, 2015.
[45] Ming Yuan and Yi Lin. Model selection and estimation in the Gaussian graphical model. Biometrika, 94(1):19–35, 2007.
[46] Xiao-Tong Yuan and Tong Zhang. Partial Gaussian graphical model estimation. IEEE Transactions on Information Theory, 60(3):1673–1687, 2014.

9

Supplement to “The Multiple Quantile GraphicalModel”
Alnur Ali
Machine Learning Department
Carnegie Mellon University
alnurali@cmu.edu

J. Zico Kolter
Computer Science Department
Carnegie Mellon University
zkolter@cs.cmu.edu

Ryan J. Tibshirani
Department of Statistics
Carnegie Mellon University
ryantibs@cmu.edu

This document contains proofs and supplementary details for the paper “The Multiple Quantile
Graphical Model”. All section, equation, and figure numbers in this supplementary document arepreceded by the letter S (all numbering without an S refers to the main paper).

S.1

Proof of Lemma 4.1

If the conditional quantiles satisfy QU |V,W (↵) = QU |W (↵) for all ↵ 2 [0, 1], then the conditional
CDF must obey the same property, i.e., FU |V,W (t) = FU |W (t) for all t in the support of U . Thisis simply because any CDF may be expressed in terms of its corresponding quantile function (i.e.,
inverse CDF), as in
FU |V,W (t) = sup{↵ 2 [0, 1] : QU |V,W (↵)  t},

and the right-hand side does not depend on V , so neither can the left-hand side. But this precisely implies that the distribution of U |V, W equals that of U |W , i.e., U and V are conditionally independent
given W . We note that the converse of the statement in the lemma is true as well, by just reversing all
the arguments here.

S.2

Proof of Lemma 4.2

This result can be seen as a generalization of Theorem 3 in [5].
First, we define an iteration of Gibbs sampling to be a single pass through all the variables (without aloss of generality, we take this order to be y1 , . . . , yd ). Now, consider a particular iteration of Gibbs
sampling; let ỹ1 , . . . , ỹd be the values assigned to the variables on the previous iteration. Then thetransition kernel for our Gibbs sampler is given byPr(y1 , . . . , yd |ỹ1 , . . . , ỹd ) = Pr(yd |yd 1 , . . . , y1 , ỹ1 , . . . , ỹd ) Pr(yd 1 , . . . , y1 |ỹ1 , . . . , ỹd ) (S.1)= Pr(yd |yd 1 , . . . , y1 ) Pr(yd 1 , . . . , y1 |ỹ1 , . . . , ỹd )(S.2)= Pr(yd |yd 1 , . . . , y1 ) Pr(yd 1 |yd 2 , . . . y1 , ỹd ) · · · Pr(y1 |ỹ2 , . . . , ỹd ),(S.3)
where (S.1) follows by the definition of conditional probability, (S.2) by conditional independence,and (S.3) by repeated applications of these tools. Since each conditional distribution is assumed
to be (strictly) positive, we have that the transition kernel is also positive, which in turn implies [2,page 544] that the induced Markov chain is ergodic with a unique stationary distribution that can bereached from any initial point.

S.3

Statement and discussion of regularity conditions for Theorem 4.3

For each k = 1, . .P
. , r, ` = 1, . . . , r, let us define the “effective” (independent) error terms(i)(i) T ⇤✏`ki = ykb⇤`k
j6=k (yj ) ✓`kj , over i = 1, . . . , n. Denote by F✏`k the conditional CDF of
29th Conference on Neural Information Processing Systems (NIPS 2016), Barcelona, Spain.

(i)✏`ki |y¬k , i = 1, . . . , n, which by construction satisfies F✏`k (0) = ↵` . Also define the underlying
support⇤S`k = j 2 {1, . . . , d} : ✓`kj
6= 0 .
Here we take a moment to explain a somewhat subtle indexing issue with the columns of the featurematrix 2 Rn⇥dm . For a single fixed index j = 1, . . . , d, we will extract an appropriate block ofcolumns of 2 Rn⇥dm , corresponding to the basis expansion of variable j, by writing j . Moreprecisely, we use j to denote the block of m columns[(j 1)m+1 ,(j 1)m+2 , . . . ,

(S.4)

jm ].

We do this because it simplifies notation considerably. (Occasionally, to be transparent, we will use
the more exhaustive notation on the right-hand side in (S.4), but this is to be treated as an exception,and the default is to use the concise notation as in j .) The same rule will be used for subsets ofindices among 1, . . . , d, so that S`k denotes the appropriate block of m|S`k | columns corresponding
to the basis expansions of the variables in S`k .
For all k = 1, . . . , d, ` = 1, . . . , r, we will assume the following regularity conditions.
c
A1. Groupwise irrepresentability: for j 2 S`k, we require that k Tj S`k kF < 1 /(6f✏`k (0) ),⇤where S`k = {j 2 {1, . . . , dm} : ✓`kj
6= 0}, f✏`k is the density of F✏`k , and > 0 is aquantity prescribed by Lemma S.5.

A2. Distributional smoothness: we assume that |F✏`k (x)|x|  C2 , where C1 , C2 > 0 are constants.

xf✏`k (0)|  C1 x2 for allF✏`k (0)A3. Correlation restriction: we assume that C3  (f✏`k (0)/n) min ( TS`k S`k )  C4 for constants C3 , C4 > 0, where min (A) denotes the minimum eigenvalue of A.
A4. Basis and support size restrictions: we assume that m = O(n1/9 ) and s = O(n1/21 ),where s = |S`k |. We also assume, with probability tending to one, that max = ⌦(1) and1/21/ log1/2 n), where we write max to denote the maximum absolute entry
max = o(nof the basis matrix .
Next, we provide some intuition for these conditions.
c
Condition A1. Fix some j 2 S`k. For notational convenience, we let

A=

T
j

S`k 2 Rm⇥sm.

Observe that each entry of A can be expressed as
Aip = n⇢i,p k(j 1)m+i k2 k

(S.5)

p k2 ,

for i = 1, . . . , m, p denoting an index into the basis expansion of the columns Sp`k , and ⇢i,p denoting
the sample correlation coefficient for the columns i and p . Since kAp kF  mkAp k1 , we havethatmax ⇢i,p <i,p1
p6n2 f✏`k (0) m

is sufficient for condition A1; here, we have also used the column scaling assumption k

p k2 

pn.

So, roughly speaking, bounded correlation between each pair of columns in the submatrices j
and S`k is enough for condition A1 to hold; note that this is trivially satisfied when Ti p = 0,
for i = 1, . . . , m, and p as defined above. Condition A1 is therefore similar to, e.g., the mutual
incoherence condition of [7] for the lasso, which is given by
T
Sc

S

T
S

1
S

1

1

˜ () maxc
j2S

T
S

1
S

T
S

j

1

1

˜,

where again S extracts the appropriate block of columns of , k · k1 here denotes the `1 operator
norm (maximum `1 norm of a row), k · k1 here denotes the elementwise `1 norm, and ˜ 2 (0, 1] is aconstant. This condition can be seen as requiring bounded correlation between each column in thesubmatrix S c and all columns in the submatrix S .
2

Condition A2. This condition is similar to requiring that f✏`k (x) be Lipschitz, over some x in aneighborhood of 0. We can show that the Laplace distribution, e.g., satisfies this condition.
The density and distribution functions for the Laplace distribution with location zero and unit scale
are given byf✏`k (x) = (1/2) exp( |x|)and⇢
1 (1/2) exp( x) if x 0F✏`k (x) =(1/2) exp(x)if x < 0,respectively.

Now, suppose 0  x  C2 . Then we can express condition A2 as|f✏`k (x)f✏`k (0)xf✏`k (0)|  C1 x2 ()2C1 x2  exp( x) + x

1  2C1 x2 .

For the first inequality, since 1 x  exp( x), it is sufficient to check that 0  C1 x2 , which is true
for C1 > 0 and all x. For the second inequality, by differentiating and again using 1 x  exp( x),we have that the function
2C1 x2 exp( x) x + 1(S.6)
is nondecreasing in x 0; thus, it is sufficient to check that this function is nonnegative for x = 0,
which is true.
Now, suppose|f✏`k (x)C2  x < 0. Then we can express condition A2 asf✏`k (0)xf✏`k (0)|  C1 x2 ()2C1 x2  exp(x)x1  2C1 x2 .

By symmetry with the preceding case, the first inequality here holds. The second inequality here also
holds, since exp(x) 2C1 x2 x 1 is continuous and increasing in x < 0; taking the limit as x " 0gives that this function is nonpositive as required.
Condition A3. This condition is a generalization of the minimum eigenvalue condition of [7], i.e.,
cmin  min (1/n) TS S , for some constant cmin > 0, and where we write S to extract theappropriate block of columns of .
Condition A4. This condition allows the number of basis functions m in the expansion to grow with n,
at a polynomial rate (with fractional exponent). This is roughly in line with standard nonparametric
regression; e.g., when estimating a continuous differentiable function via a spline expansion, one
typically takes the number of basis functions m to scale as n1/3 [4]. The condition also restricts,
for any given variable, the number of variables s that contribute to its neighborhood model to bepolynomial in n (with a smaller fractional exponent).
Finally, the condition assumes that the entries of the basis matrix (i.e., the matrix of transformed
variables) to be at least of constant order, and at most of polynomial order (with small fractional
exponent), with n. We note that this implicitly places a restriction on the tails of distribution governing(i)the data yj , i = 1, . . . , n, j = 1, . . . , d. However, the restriction is not a strong one, because it
allows the maximum to grow polynomially large with n (whereas a logarithmic growth would be
expected, e.g., for normal data). Furthermore, it is possible to trade off the restrictions on m, s, max ,and d (presented in the statement of the theorem), making each of these restrictions more or lessstringent, if required.

S.4

Proof of Theorem 4.3

The general strategy that we use here for support recovery is inspired by that in [3], for `1 -penalized
quantile regression.
Fix some k = 1, . . . , d and ` = 1, . . . , r. We consider the conditional distribution yk |y¬k , whose↵` -quantile is assumed to satisfy (3). Hence, to be perfectly clear, all expectations and probability(i)statements in what follows are to be interpreted with respect to the observations yk , i = 1, . . . , n(i)conditional on yj , i = 1, . . . , n, for j 6= k (and thus we can treat the feature matrix as fixed3

throughout). In the setting assumed by the theorem, the conditional quantile model in (3) is, more
explicitly,dX⇤
Qyk |y¬k (↵` ) = b⇤`k +(✓`kj)T j (yj ),j6=k⇤for some unknown parameters b⇤`k and ✓`kj
, j = 1, . . . , d. For simplicity, in this proof, we will drop
the intercept term completely both from the model (denoted b⇤`k ) and the optimization problem in (4)(here denoted b`k ) that defines the estimator in question. Including the intercept is not at all difficult,and it just requires some extra bookkeeping at various places. Recall that we define⇤S`k = j 2 {1, . . . , d} : ✓`kj
6= 0 ,

and analogously defineŜ`k = j 2 {1, . . . , d} : ✓ˆ`kj 6= 0 ,ˆˆˆwhere ✓`k = (✓`k1 , . . . , ✓`kd ) 2 Rdm is the solution in (5).

We will show that, with probability at least 1/(dr), it holds that S`k = Ŝ`k . A union bound (overall choices k = 1, . . . , d and ` = 1, . . . , r) will then tell us that E ⇤ = Ê with probability at least1, completing the proof.
To certify that S`k = Ŝ`k , we will show that the unique solution in (5) is given by✓ˆ`k(S`k ) = ✓˜`k(S`k ) ,c ) = 0,✓ˆ`k(S`k(S.7)

where ✓˜`k(S`k ) solves the “restricted” optimization problem:⇣⌘
Xminimize ↵` Yk✓k✓`kj k2 +S`k `k(S`k ) + 1✓`k(S`k )

j2S`k

2

2

k✓`k(S`k ) k22 .

(S.8)

Now, to prove that ✓ˆ`k as defined above in (S.7) indeed the solution in (5), we need to check that it
satisfies the KKT conditions for (5), namely
⇣
⌘
T
˜
˜
(S.9)
S`k ✓`k(S`k )2 ✓`k(S`k ) = 1 u`k(S`k ) ,S`k v` Yk
⇣
⌘
T
˜
c ),
(S.10)
c v ` YkS`k ✓`k(S`k ) = 1 u`k(S`kS`kn˜where v` (YkS`k ✓`k(S`k ) ) 2 R is a subgradient of
h ⇣
⌘i
⇣
(i)
˜
v` Yk
= ↵ ` I yk
S`k ✓`k(S`k )i˜↵` (·) at YkS`k ✓`k(S`k ) , i.e.,˜i(S`k ) ✓`k(S`k )⌘,i = 1, . . . , n

where I (·) is the indicator function of the nonpositive real line, and where each u`kj 2 Rm is a
subgradient of k · k2 at ✓˜`kj , i.e.,⇢{✓˜`kj /k✓˜`kj k2 }if ✓`kj 6= 0u`kj 2{x 2 Rm : kxk2  1} if ✓`kj = 0,

for j = 1, . . . , d. Note that, since ✓˜`k(S`k ) is optimal for the restricted problem (S.8), we know that
there exists a collection of subgradients u`k(S`k ) to satisfy (S.9), from the KKT conditions for (S.8)
itself.
˜
It remains to satisfy (S.10), and for this, we can use u`kj = Tj v` (YkS`k ✓`k(S`k ) ) as a validcchoice of subgradient, for each j 2 S`k , provided that
⇣
⌘
T
c
˜
< 1 , for j 2 S`k.
(S.11)
S`k ✓`k(S`k )j v` Yk
Define zj (#) =Tj v` (Yk2cS`k #), for j 2 S`k , and define a ball⇤B ⇤ = # 2 Rsm : k# ✓`k(Sk `k ) 2,where we write s = |S`k |. To show (S.11), then, it suffices to show that✓˜`k(S`k ) 2 B ⇤ ,
|
{z
}
E1

and

max sup kzj (#)k2 < 1 .
|
{z
}
c
j2S`k
#2B ⇤

E2

4

(S.12)

In Lemma S.5.1, given in Section S.5, it is shown that the event E1 defined above occurs withprobability at least 1/(2dr), with a choice of radius!
rps log n
1s m
=C
+
,
n
n
for a constant C > 0. Below we show that E2 occurs with probability at least 1
For j = 1, . . . , d, let us expandh ⇣zj (#) = Tj v` (✏`k ) + Tj E v` Yk| {z } |
j
1

|

✓ ⇣
T
j v` Yk

S`k #

{z

j
2

S`k #)⌘v` (✏`kiv` (✏`k ) +
}

⌘

{z

j
3

h
E v` (YkS`k #)/(2dr), as well.

v` (✏`k )i◆,}(S.13)

where ✏`k = (✏`k1 , . . . , ✏`kn ) 2 Rn is a vector of the effective error terms, which recall, is defined by⇤✏`k = Yk✓`k
. Therefore, to show that the event E2 in (S.12) holds, we can show that for eachp = 1, 2, 3,
max sup k

c
j2S`k
#2B ⇤

j
p k2 <

1

3.

Further, to show that E2 holds with probability at least 1/(2dr), we can show that the above holds
for p = 1, 3 each with probability at least 1/(4dr), as the statement for p = 2 is deterministic.
We now bound the terms j1 , j2 , j3 one by one.
Bounding kjc
1 k2 . Fix j 2 S`k , and writeTj v` (✏`k ) =✓Xni,(j 1)m+1 v` (✏`ki ), . . . ,

i=1

n
X

i,jm v` (✏`ki )

i=1

◆

,

where, as a reminder that the above quantity is a vector, we have returned momentarily to the moreexhaustive notation for indexing the columns of , as in the right-hand side of (S.4).
Straightforward calculations reveal that, for each i = 1, . . . , n, and p = 1, . . . , m,Ei,(j 1)m+p v` (✏`ki ) = 0,and|Hence,⇣Pr kTj v` (✏`ki )k2p⌘mt  Pri,(j 1)m+p | 

✓ X
n

m
X

i,(j 1)m+p v` (✏`ki )  |i,(j 1)m+p v` (✏`ki )t, some p = 1, . . . , m

i=1

2 exp

p=1

 2m exp

✓

✓

Pn

2 i=1
◆
t2
.
2n

i,(j 1)m+p |.

t2
2i,(j 1)m+p◆◆p
Above, the first inequality used the simple fact that kxk2  mkxk1 for x 2 Rm ; the second used
Hoeffding’s bound p
and the union bound; and the thirdp used our assumption that the columns ofhave norm at most n. Therefore, taking t = 1 /(3 m), we see that, by the above and the union
bound,✓◆✓◆
Pr

By choosing 1 = C 0
question is at least 1

max
k
c

j2S`k

j
1 k2 <

1

1

3

2dm exp

2
1

18mn

.

p
18mn log(8d2 mr/ ) for a constant C 0 > 0, we see that the probability in/(4dr), as desired.
5

Bounding k j2 k2 . Recall that F✏`k (·) is used to denote the CDF of the effective error distribution,and f✏`k (·) is used for its its density. By construction, F✏`k (0) = ↵` . Direct calculation, using thedefinition of v` (·), shows that, for any # 2 B ⇤ , and each i = 1, . . . , n,
h
⇣
⌘i
⇣
⌘
⇤
E v` (✏`k ) v` Yk= F✏`k S`k # ✓`k(SF✏`k (0),
S`k #
`k )

where we apply F✏`k componentwise, and so
h
⇣
⌘i
T
E
v
(✏
)
v
Y
#
= f✏`k (0)
`
`k
`
k
S
j
`k
with

h

⇤
✓`k(S
+
`k )

F✏`k (0)f✏`k (0)i(S`k ) #

T
j

#

⇤
✓`k(S

`k ) 2S`kjm
4 2 R being the appropriate remainder term, i.e.,

j
4

i

t

=

n
X

it

i=1

for t = j(m

h

F✏`k

⇣

⇤
✓`k(S
`k )

i(S`k ) #

⌘

j
4

#

T
j

i
⇤
✓`k(S
,
`k )

1) + 1, . . . , jm.

Now, we have thatf✏`k (0)Tjwhere we have used k#

⇤
✓`k(S
`k )

#

S`k

⇤
✓`k(S
k2 
`k )

2

 f✏`k (0)

S`k F

1

6

,

and the groupwise irrepresentability condition in A1.

We also have the following two facts, which we will use momentarily:32= o( 1 )max nsps(S.14)(S.15)! 0.

max

Note that (S.14) can be obtained as follows. Since (1/2)(x + y)2  x2 + y 2 for x, y 2 R, we canplug in!
rps log n
1s m
=C
+
,
n
n
and check that both terms on the right-hand side of
✓ 2 2
◆
3
s log n
max ns
1s m
+
=
n2
n
1

3
3
max s 1 mn

3
2
max s log n+1tend to zero. For the first term on the right-hand side, it is enough to show that
3
6
6 3
2
max s m log(d mr)(log n)/n ! 0,q
where we have plugged in 1 = C 0 mn log(d2 mr/ ) log3 n. Using the assumptions in condition
A4, we get that log(d2 mr) = O(log d + log m) = O(n2/21 ), and furthermore that✓ 1/3 2/21 6/21 6/21 ◆
n ·n
·n
·n
log3 n
3
6
6 3
2
! 0,
max s m log(d mr)(log n)/n = onlog3 nas required. A similar calculation shows that the second term on the right-hand side also tends tozero, i.e., 3max s2 (log n)/ 1 ! 0, which establishes (S.14). Lastly, (S.15) follows since its left-handside is dominated by the left-hand side of (S.14).
So, we now compute
k

p
j
m max
4 k2 
t

n
X
i=1

it

h

F✏`k

⇣

i(S`k )

#

⇤
✓`k(S
`k )

⌘
F✏`k (0)

 C1

max

p

m

n ⇣
X

i(S`k )

#

⇤
✓`k(S
`k )

i=1

6

⌘2

f✏`k (0)i(S`k ) #

⇤
✓`k(S
`k )

i

 C1

max

p

m

n
X
i=1

k

2
i(S`k ) k2 k#

⇤
✓`k(S
k2
`k ) 2

p
 C1 3max mns 2= o( 1 ).
p
Here the first inequality follows from the fact that kxk2  mkxk1 for x 2 Rm , and the triangle
inequality; the second follows from the distributional smoothness condition in A2, which is applicablesince (S.15) holds; the third uses Cauchy-Schwarz; the fourth uses our column norm assumption on⇤, and k# ✓`k(Sk2  ; the last uses (S.14). As k j4 k2 = o( 1 ), it will certainly be strictly less`k )than 1 /6 for n large enough. We have hence shown, noting that none of our above arguments havedepended on the particular choice of j = 1, . . . , d or # 2 B ⇤ ,max sup k

c
j2S`k
#2B ⇤

j
2 k2 <

1

3.

Bounding k j3 k2 . For this part, we can use the end of the proof of Lemma 2 in [3], which uses classic
entropy-based techniques to establish a bound very similar to that which we are seeking. By carefully
looking at the conditions required for this lemma, we see that under the distributional smoothness
condition in A2, condition A3, and alsopn log(dm) = o( 1 )
n

(1 +

max

2

= o( 1 )2
3/2) log2 n = o( 21 /n),max s

all following directly from condition A4 by calculations similar to the ones we used when bounding
k j2 k, we have
✓
◆
✓
◆
1
1
j
j
p
Pr max
sup k 3 k2
 Pr max
sup k 3 k1
;
c
c
j2S`k
j2S`k
3
3 m
#2B ⇤
#2B ⇤
the probability on the right-hand side can be made arbitrarily small for large n, by the arguments at
the end of Lemma 2 in [3], and hence clearly smaller than the desired /(4dr) level.
Putting it together. Returning to the logic in (S.11), (S.12), (S.13), we have shown that the subgradientcondition in (S.11) holds with probability at least 1 ( /(2dr) + /(4dr) + /(4dr)) = 1/(dr).
Taking a union bound over k = 1, . . . , d and ` = 1, . . . , r, which were considered fixed at the start ofour analysis, gives the result stated in the theorem.

S.5

Statement and proof of Lemma S.5.1

We show that with probability at least 1/(2dr), it holds that ✓˜`k(S`k ) 2 B ⇤ , where ✓˜`k(S`k ) is thesolution to the restricted problem (S.8), for some fixed k = 1, . . . , d and ` = 1, . . . , r, and B ⇤ is aball defined in the proof of Theorem 4.3 in Section S.4. This fact is used a few times in the proof of
Theorem 4.3.
Lemma S.5.1. Fix some k = 1, . . . , d and ` = 1, . . . , r. Define the ballB ⇤ = {# 2 Rsm : k#

⇤
✓`k(S
k  }
`k ) 2⇤
centered at the underlying coefficients ✓`k(Swith radius
`k )
!
r
p
s log n
1s m
=C
+
,
n
n

for some constant C > 0. Then, with probability at least 1
where ✓˜`k(S`k ) is the solution to the restricted problem (S.8).
7/(2dr), it holds that ✓˜`k(S`k ) 2 B ⇤ ,Proof. We will follow the strategy for the proof of Theorem 1 in [3] closely. We begin by considering
the ball⇤B = {# 2 Rsm : k# ✓`k(Sk  R}`k ) 2⇤
with center ✓`k(Sand radius R. We also introduce some useful notational shorthand, and write the`k )quantile loss term in the restricted problem (S.8) asL`k (#) =↵` (YkS`k #) .

Below, we show that a particular function of R serves as an upper bound for the quantity⇤
E[L`k (#̃`k(S`k ) ) L`k (✓`k(S)], where the expectation here is taken over draws of the data, and`k )#̃`k(S`k ) is a particular point in B that we define in a moment. This in turn implies, with probability
at least 1/(2dr), that ✓˜`k(S`k ) 2 B ⇤ , as claimed.
First, we define #̃`k(S`k ) more precisely: it is a point on the line segment between the solution to the⇤restricted problem ✓˜`k(S`k ) and the underlying coefficients ✓`k(S, i.e.,`k )#̃`k(S`k ) = ✓˜`k(S`k ) + (1for a particular choice⇤)✓`k(S
,
`k )

R
,
⇤
R + k✓˜`k(S`k ) ✓`k(S
k2
`k )
˜
/ B , as we establish next. Observe that we) 2 B even if ✓`k(S ) 2=

which guarantees that #̃`k(S`kalways havek✓˜`k(S`k )
() R
()

`k

⇤
✓`k(S
k  R + k✓˜`k(S`k )`k ) 2⇤k✓˜`k(S`k ) ✓`k(S
k2
`k )
R
R + k✓˜`k(S ) ✓⇤k2k✓˜`k(S`k )() k ✓˜`k(S`k )() k#̃`k(S`k )

`k

⇤
✓`k(S
k
`k ) 2`k(S`k )

⇤
✓`k(S
k R
`k ) 2⇤⇤✓`k(S+ ✓`k(S
`k )
`k )

⇤
✓`k(S
k  R,`k ) 2

⇤
✓`k(S
k R
`k ) 2

as claimed. The second line here follows by rearranging and multiplying through by R; the third byusing the definition of above; the fourth by adding and subtracting the underlying coefficients; andthe fifth by using the definition of #̃`k(S`k ) .
Now, the beginning of the proof of Theorem 1 in [3] establishes, for any #̃`k(S`k ) 2 B, for someconstant C5 > 0, and using condition A3, thathi⇤⇤
E L`k (#̃`k(S`k ) ) L`k (✓`k(S)
C5 nk#̃`k(S`k ) ✓`k(Sk2 ,(S.16)`k )`k ) 2and so, by direct calculation, since⇤⇤k✓˜`k(S`k ) ✓`k(Sk  R () k✓˜`k(S`k ) ✓`k(Sk  R/2,`k ) 2`k ) 2(S.17)⇤
it suffices to obtain a suitable upper bound for E[L`k (#̃`k(S`k ) ) L`k (✓`k(S
)],
in
order
to
get
the
`k )
result in the statement of the lemma. To this end, we introduce one more piece of shorthand, anddenote the objective for the restricted problem (S.8) as J`k (#).
⇤
k#̃`k(S`k ) ✓`k(Sk  R ()`k ) 2

We proceed with the following chain of (in)equalities:hi⇤
E L`k (#̃`k(S`k ) ) L`k (✓`k(S
)
)
`k
h
i
⇤
= E L`k (#̃`k(S`k ) ) L`k (✓`k(S`k ) ) + J`k (#̃`k(S`k ) )8J`k (#̃`k(S`k ) ) +⇤J`k (✓`k(S
)
`k )

⇤
J`k (✓`k(S)`k )(S.18)

⇤= L`k (✓`k(S
)
`k )
|

J`k (#̃`k(S`k ) )⇤E L`k (✓`k(S)`k )⇤(✓`k(S{z`k )L`k (#̃`k(S`k ) ) + E L`k (#̃`k(S`k ) ) +},#̃`k(S`k ) )⇤J`k (✓`k(S
)+
`k )

1

X

j2S`k



⇤
(✓`k(S, #̃`k(S`k ) ) + J`k (#̃`k(S`k ) )
`k )



⇤
(✓`k(S, #̃`k(S`k ) ) + J`k (#̃`k(S`k ) )
`k )




⇤
k✓`kj
k2

X

1

j2S`k

k#̃`kj k2⇤( 2 /2)k#̃`k(S`k ) k22 + ( 2 /2)k✓`k(Sk2`k ) 2(S.19)X⇤⇤J`k (✓`k(S`k ) ) + 1
k✓`kj #̃`kj k2 + o(1)
j2S`k

p

⇤
J`k (✓`k(S
)+
`k )

⇤
(✓`k(S, #̃`k(S`k ) ) + 2 1 sR m
`k )

sup
#̃`k(S`k ) 2B|⇤(✓`k(S, #̃`k(S`k ) )| + 2 1 sR
`k )

p

1 sR

p

(S.20)(S.21)

m + o(1)(S.22)(S.23)

m.

Here, (S.18) follows by adding and subtracting like terms, and (S.19) by rearranging (S.18). In (S.20)
we use the triangle inequality and the following argument to show that the terms involving 2 are o(1).
41/42 ⇤
Under the assumption/✓max ), combined with the restriction that s = o(n1/21 ),2 = o(np that⇤we have 2 = o(n/( s✓max )). Therefore, under our choice of R = 1/n (as specified below), wehavep ⇤2 s✓max R ! 0.
This in turn is used to argue that⇤( 2 /2)k#̃`k(S`k ) k22 + ( 2 /2)k✓`k(Sk2 = ( 2 /2)k#̃`k(S`k )`k ) 2

 ( 2 /2)R2

 ( 2 /2)R2 +2

 ( 2 /2)R + ( 2 /2)R2 +

⇤
✓`k(S
k2
`k ) 22
2 k#̃`k(S`k ) k2 +T⇤
2 #̃`k(S`k ) ✓`k(S`k )2 k#̃`k(S`k ) k2 (k#̃`k(S`k ) k2⇤k✓`k(Sk )`k ) 2

2 k#̃`k(S`k ) k2 R⇤
2 k✓`k(S`k ) k2 R
2

p

⇤
s✓max
R ! 0.

⇤In the second to last line, we have applied k#̃`k(S`k ) k2  k✓`k(Sk2 , as, outside of this case, the`k )⇤2term in question ( 2 /2)k#̃`k(S`k ) k22 + ( 2 /2)k✓`k(S
k
would
be
negative, anyway.
`k ) 2⇤⇤
Continuing on, (S.21) holds because k✓`k(S#̃`k(S`k ) k2  R implies k✓`kj#̃`kj k2  R. Fi`k )nally, (S.22) follows because of the following argument. Since J`k is convex, we can use the definitionof #̃`k(S`k ) and getJ`k (#̃`k(S`k ) )  J`k (✓˜`k(S`k ) )+(1⇤⇤⇤⇤)J`k (✓`k(S) = J`k (✓`k(S)+ (J`k (✓˜`k(S) J`k (✓`k(S
));
`k )
`k )
`k )
`k )

notice that the last term here is nonpositive, since ✓˜`k(S`k ) is the solution to the restricted problem(S.8), and thus we have that⇤J`k (#̃`k(S`k ) )  J`k (✓`k(S),`k )which lets us move from (S.21) to (S.22).

Lemma 1 in [3] states, with probability at least 1
some constant, thatsup#̃`k(S`k ) 2B

, where⇤, #̃`k(S`k ) )|  6R
| (✓`k(S`k )= exp( C6 s log n) and C6 > 0 isp

sn log n,

so from (S.23), with probability at least 1, we see that
h
i
p
p
⇤
E L`k (#̃`k(S`k ) ) L`k (✓`k(S)
 6R sn log n + 2 1 sR m`k )9and, using (S.16), thatnk#̃`k(S`k )for some constant C 0 > 0.

⇣ p
⇤
2
0
✓`k(S
k

C
R sn log n +`k ) 2

p ⌘
sR
m ,
1

Plugging in R = 1/n, dividing through by n, and using the fact that the square root function is
subadditive, we get, with probability at least 1, that✓◆(s log n)1/4( 1 s)1/2 m1/4⇤
0
k#̃`k(S`k ) ✓`k(S`k ) k2  C
+
n
n3/4
r
p !
s log n
1s m
 C0
+
.
n
n
Finally, we complete the proof by applying (S.17), in order to get thatk✓˜`k(S`k )where we have defined

=C

⇤
✓`k(S
k  ,
`k ) 2

1s

p

m

n

+

r

and C > 0 is some constant, with probability at least 1

S.6

s log nn!

,/(2dr), for large enough n.

Proof of Lemma 5.1

The prox operator prox A (A) is separable in the entries of its minimizer X, so we can focus on
minimizing over Xij the expression
max{↵j Xij , (↵j1)Xij } + (1/(2 )) (Xij= ↵j max{0, Xij } + (1

Aij )

2

↵j ) max{0, Xij } + (1/(2 )) (Xij2

Aij ) .

(S.24)

Suppose Xij > 0. Then differentiating (S.24) gives Xij = Aij↵j and the sufficient condition
Aij > ↵j . Similarly, assuming Xij < 0 gives Xij = Aij + (1 ↵j ) when Aij < (↵j 1).
Otherwise, we can take Xij = 0. Putting these cases together gives the result.

S.7

ADMM for the MQGM

A complete description of our ADMM-based algorithm for fitting the MQGM to data is given in
Algorithm 1.

S.8

Additional details on Gibbs sampling

In the MQGM, there is no analytic solution for parameters like the mean, median, or quantiles ofthese marginal and conditional distributions, but the pseudolikelihood approximation makes for avery efficient Gibbs sampling procedure, which we highlight in this section. As it is relevant tothe computational aspects of the approach, in this subsection we will make explicit the conditionalrandom field, where yk depends on both y¬k and fixed input features x.
First, note that since we are representing the distribution of yk |y¬k , x via its inverse CDF, to sample
from from this conditional distribution we can simply generate a random ↵ ⇠ Uniform(0, 1). We
then computexQ̂yk |y¬k (↵` ) = (y)T ✓`k + xT ✓`kxQ̂yk |y¬k (↵`+1 ) = (y)T ✓(`+1)k + xT ✓(`+1)kfor some pair ↵`  ↵  ↵`+1 and set yk to be a linear interpolation of the two values,
10

Algorithm 1 ADMM for the MQGM
Input: observations y (1) , . . . , y (n) 2 Rd , feature matrix
1, 2 > 0
ˆ = (✓ˆ`kj , b̂`k )Output: fitted coefficients ⇥for k = 1, . . . , d (in parallel, if possible) do
initialize ⇥k , Bk , V, W, Z, UV , UW , UZ
repeatupdate ⇥k using (13)update Bk using (13)update V using (10)update W using (11)update Z using (12) and Lemma 5.1
update UV , UW , UZ :2 Rn⇥dm , quantile levels A, constants

UVUW

UV + (1BkT + k ⇥UW + (⇥k W )UZUZ + (Yk 1T

V)

1BkT

k⇥

Z)

until converged
end foryk

Q̂yk |y¬k (↵` ) +⇣Q̂yk |y¬k (↵`+1 )⌘Q̂yk |y¬k (↵` ) (↵

↵`+1

↵`

↵` )

.

This highlights the desirability of having a range of nonuniformly spaced ↵ terms that reach values
close to zero and one as otherwise we may not be able to find a pair of ↵’s that lower and upper boundour random sample ↵. However, in the case that we model a sufficient quantity of ↵, a reasonableapproximation (albeit one that will not sample from the extreme tails) is also simply to pick a random↵` 2 A and use just the corresponding column ✓`k to generate the random sample.

Computationally, there are a few simple but key points involved in making the sampling efficient.
First, when sampling from a conditional distribution, we can precompute xT ⇥xk for each k, and use
these terms as a constant offset. Second, we maintain a “running” feature vector (y) 2 Rdm , i.e.,the concatenation of features corresponding to each coordinate (yk ). Each time we sample a newcoordinate yk , we generate just the new features in the (yk ) block, leaving the remaining features
untouched. Finally, since the ⇥k terms are sparse, the inner product (y)T ✓`k will only contain a
few nonzeros terms in the sum, and will be computed more efficiently if the ⇥k are stored as a sparse
matrices.

S.9

Additional numerical results for the ring data

S.9.1

Conditional independencies recovered by the nonparanormal skeptic, TIGER, andLaplace

We present the conditional independencies recovered by the nonparanormal skeptic, TIGER, andLaplace on the ring data in Figure S.1; results for the remaining methods are presented in Section 6.1.
S.9.2

Evaluation of fitted conditional CDFs

Here, we elaborate on the evaluation of the conditional CDFs given by the MQGM, MB, GLasso,
SpaceJam, TIGER, and Laplace, when run on the ring data. (We omit the nonparanormal skeptic fromour evaluation as it is not clear how to sample from its conditionals, due to the nature of a particulartransformation that it uses.)
For each of these methods, we essentially averaged the total variation distances and KolmogorovSmirnoff statistics between the fitted and true conditional CDFs across all variables, and then reported11

λ1 =0.00781

λ1 =0.01562

λ1 =0.03125

λ1 =0.06250

λ1 =0.12500

Truth

λ1 =0.03125

λ1 =0.06250

λ1 =0.12500

λ1 =0.25000

λ1 =0.50000

Truth

λ1 =4.00000

λ1 =8.00000

λ1 =16.00000

λ1 =32.00000

λ1 =64.00000

Laplace

TIGE5

1onpara.

7ruth

Figure S.1: Conditional independencies recovered by the nonparanormal skeptic, TIGER, and Laplace on thering data; black means conditional dependence.
Table S.1: Total variation (TV) distance and Kolmogorov-Smirnoff (KS) statistic values for the MQGM, MB,
GLasso, SpaceJam, TIGER, and Laplace on the ring data; lower is better, best in bold.

MQGM
MB
GLassoSpaceJam
TIGER
LaplaceTV
20.87392.29892.47991.56888.284
127.406

KS
0.760
1.856
1.768
1.697
1.450
1.768

the best values obtained across a range of tuning parameters (more details below). We present theresults in Table S.1; we see that the MQGM outperforms all its competitors, in both metrics.
We now describe the evaluation in more detail; for simplicity, we describe everything that follows interms of the conditional CDF y1 |y2 only, with everything being extended in the obvious way to otherconditionals.
First, we carried out the following steps in order to compute the true (empirical) conditional CDF.
1. We drew n = 400 samples from the ring distribution, by following the procedure describedin Section 6.1; these observations are plotted across the top row of Figure S.2.
2. We then partitioned the y2 samples into five equally-sized bins, and computed the true
empirical conditional CDF of y1 given each bin of y2 values.
Next, we carried out the following steps in order to compute the estimated (empirical) conditional
CDFs, for each method.
3. We fitted each method to the samples obtained in step (1) above.
4. Then, for each method, we drew a sample of y1 given each y2 sample, using the method’s
conditional distribution; these conditionals are plotted across the second through fifth rowsof Figure S.2 (for representative values of 1 ).
Operationally, we drew samples from each method’s conditionals in the following ways.
• MQGM: we used the Gibbs sampler described in Section S.8.
(i)
• MB: we drew y1 ⇠ N (✓ˆT y , ˆ 2 ), where ✓ˆ1 is the fitted lasso regression coefficient
1

2

1|2

(i)

of y1 on y2 ; y2 for i = 1, . . . , n is the ith observation of y2 obtained in step (1) above;2and ˆ1|2= var(Y1 Y2 ✓ˆ1 ) denotes the sample variance of the underlying error term
(1)(n)
Y1 Y2 ✓ˆ1 with Yi = (yi , . . . , yi ) 2 Rn collecting all observations along variable
i.
12

(i)
2
• SpaceJam: we drew y1 ⇠ N (✓ˆ1T (y2 ), ˆ1|2), where is a suitable basis function,2ˆand ✓1 as well as ˆ1|2 are defined in ways analogous to the neighborhood selection
setup.
2• GLasso: we drew y1 ⇠ N (µ̂1|2 , ˆ1|2), whereˆ 12 ⌃
ˆ 1 (y (i)µ̂1|2 = µ̂1 + ⌃
22
2
1ˆ
2
ˆ
ˆ
ˆ
ˆ = ⌃11 ⌃12 ⌃ ⌃21

µ̂2 )

22

1|2

ˆ denoting the estimate of the covariancewith µ̂i denoting the sample mean of Yi , and ⌃matrix given by GLasso (subscripts select blocks of this matrix).
5. Finally, we partitioned the y2 samples into five equally-sized bins (just as when computing
the true conditional CDF), and computed the estimated empirical conditional CDF of y1
given each bin of y2 values.
Having computed the true as well as estimated conditional CDFs, we measured the goodness of fit ofeach method’s conditional CDFs to the true conditional CDFs, by computing the total variation (TV)distance, i.e.,qXmethod(z (i) |⇣) ,(1/2)
F̂y1 |y2 j (z (i) |⇣) F̂ytrue
1 |y2
i=1

as well as the (scaled) Kolmogorov-Smirnoff (KS) statistic, i.e.,max

i=1,...,q

methodF̂y1 |y2 j (z (i) |⇣)F̂ytrue(z (i) |⇣) .
1 |y2(z (i) |⇣) is the true empirical conditional CDF of y1 |y2 , evaluated at y1 = z (i) and
Here, F̂ytrue
1 |y2methodgiven y2 = ⇣, and F̂y1 |y2 j (z (i) |⇣) is a particular method’s (“methodj ” above) estimated empiricalconditional CDF, evaluated at y1 = z (i) and given y2 = ⇣. For each method, we averaged these TV
and KS values across the method’s conditional CDFs. Table S.1 reports the best (across a range oftuning parameters) of these averaged TV and KS values.

S.10

Additional numerical results for modeling flu epidemics

Here, we plot samples from the marginal distributions of the percentages of flu reports at regions
one, five, and ten throughout the year, which reveals the heteroskedastic nature of the data (just as in
Section 6.2, for region six).

S.11

Sustainable energy application

We evaluate the ability of MQGM to recover the conditional independencies between several windfarms on the basis of large-scale, hourly wind power measurements; wind power is intermittent, andthus understanding the relationships between wind farms can help farm operators plan. We obtainedhourly wind power measurements from July 1, 2009 through September 14, 2010 at seven wind farms(n = 877, see [6, 8, 1] for details). The primary variables here encode the hourly wind power at a
farm over two days (i.e., 48 hours), thus d = 7 ⇥ 48 = 336. Exogenous variables were used to encode
forecasted wind power and direction as well as other historical measurements, for a total of q = 3417.
We set m = 5 and r = 20. Fitting the MQGM here hence requires solving 48 ⇥ 7 = 336 multiple
quantile regression subproblems each of dimension ((336 1) ⇥ 5 + 3417) ⇥ 20 = 101, 840. Each
subproblem took roughly 87 minutes, comparable to the algorithm of [8].
Figure S.4 presents the recovered conditional independencies; the nonzero super- and sub-diagonal
entries suggest that at any wind farm, the previous hour’s wind power (naturally) influences thenext hour’s, while the nonzero off-diagonal entries, e.g., in the (4,6) block, uncover farms that may
influence one another. [8], whose method placed fifth in a Kaggle competition, as well as [1] reportsimilar findings (see the left panels of Figures 7 and 3 in these papers, respectively).

13

y1 bin 1

1.5

y1 bin 3

1.5

truth
1.0

0.5

0.5

0.5

0.0

0.0

0.0

y1

1.0

y1

y1

truth

1.0

−0.5

−0.5

−0.5

−1.0

−1.0

−1.0

−1.0

−0.5

0.0

0.5

1.0

1.5

−0.5

0.0

0.5

1.0

1.5

04G0

1.0

1.5

0.0
−1.5

0.6

0.2

−1.0

−0.5

0.0

0.5

1.0

1.5

0.0
−1.5

y1

y1 bin 3

y1 bin 5
0B

0.6

0.5

1.0

1.5

0.0
−1.5

0.2

−1.0

−0.5

y1

0.0

0.5

1.0

1.5

0.0
−1.5

y1 bLn 1−0.5

Pr(y1 |y1 ∈ [0.790,1.294])

3r(y1 |y1 ∈ [-0.217,0.287])

Pr(y1 |y1 ∈ [-1.224,-0.721])

0.6

0.4

0.4

0.2

0.2

0.2

0.0
−1.5

0.0
−1.5

0.0
−1.5

0.5

1.0

1.5

−1.0

−0.5

0.0

0.5

1.0

1.5

−1.0

−0.5

0.0

y1

y1

y1

y1 bin 1

y1 bin 3

y1 bin 5

6Sace-aP

6Sace-aP

1.0

0.8

0.6

0.4

0.2

0.0

y1

0.5

1.0

1.5

0.0
−1.5

6SaceJaP

0.6

0.4

0.2

1.5

0.8

0.6

0.4

1.0

Pr(y1 |y1 ∈ [0.790,1.294])

0.8

0.5

1.0

3r(y1 |y1 ∈ [-0.217,0.287])

Pr(y1 |y1 ∈ [-1.224,-0.721])

1.0

1.5

0.8

0.6

0.4

1.0

GLasso

1.0

0.8

0.6

0.5

y1 bLn 5
GLasso

1.0

0.8

0.0

y1

y1 bLn 3
GLasso

−0.5

−1.0

y1

1.0

0.0

0B

0.4

0.2

0.0

1.5

0.6

0.4

0.2

1.0

0.8

0.6

0.4

0.5

1.0

3r(y1 |y1 ∈ [-0.217,0.287])

Pr(y1 |y1 ∈ [-1.224,-0.721])

1.0

0.8

−1.0

0.0

y1 bin 1

0.8

0.0
−1.5

−0.5

y1

0B

−0.5

−1.0

y1

1.0

−1.0

04G0

0.4

0.2

0.5

1.5

Pr(y1 |y1 ∈ [0.790,1.294])

0.6

0.0

1.0

0.8

0.4

0.2

0.5

1.0

3r(y1 |y1 ∈ [-0.217,0.287])

Pr(y1 |y1 ∈ [-1.224,-0.721])

1.0

0.6

−0.5

0.0

y1

0.4

−1.0

−0.5

y1 biQ 5

0.8

−0.5

−1.0

y1 biQ 3
04G0

−1.0

−1.5
−1.5

y1 biQ 1

0.8

0.0
−1.5

−1.0

y1

1.0

0.0
−1.5

−1.5
−1.5

y1

Pr(y1 |y1 ∈ [0.790,1.294])

−1.5
−1.5

y1 bin 5

1.5

truth

0.2

−1.0

−0.5

0.0

y1

0.5

1.0

1.5

0.0
−1.5

−1.0

−0.5

0.0

0.5

1.0

1.5

y1

Figure S.2: Conditional distributions for MQGM, MB, GLasso, and SpaceJam, fitted to samples from the ring
distribution (TIGER and Laplace’s conditionals both look similar to MB’s). First row: samples from the ring
distribution, where each plot highlights the samples falling into a particular shaded bin on the y2 axis. Second
through fifth rows: conditional distributions of y1 given y2 for each method, where each plot conditions on theappropriate y2 bin as highlighted in the first row. The MQGM’s conditional distributions are intuitive, appearing
bimodal for bin 3, and more peaked for bins 1 and 5. MB, GLasso, and SpaceJam’s densities appear (roughly)
Gaussian, as expected.

14

5egiRn 1

7
6
5
4
3
2
1
40

50

8

18

28

8
7
6
5
4
3
2
1
0
30

40

Week

50

8

5egiRn 99% Rf flu-like sympWRms

8

0
30

5egiRn 3

9% Rf flu-like sympWRms% Rf flu-like sympWRms

9

18

28

8
7
6
5
4
3
2
1
0
30

Week

40

50

8

18

28

Week

Figure S.3: Samples from the fitted marginal distributions of the weekly flu incidence rates at several regions of
the U.S.; samples at larger quantile levels shaded lighter, median in darker blue.

0

0

48

96

144

192

240

288

336

48
96
144
192
240
288
336

Figure S.4: Conditional independencies recovered by the MQGM on the wind farms data; each block corresponds to a wind farm, and black indicates dependence.

15

References[1] Alnur Ali, Kshitij Khare, Sang-Yun Oh, and Bala Rajaratnam. Generalized pseudolikelihood
methods for inverse covariance estimation. Technical report, 2016. Available at http://arxiv.
org/pdf/1606.00033.pdf.
[2] Christopher Bishop. Pattern Recognition and Machine Learning. Springer, 2006.
[3] Jianqing Fan, Yingying Fan, and Emre Barut. Adaptive robust variable selection. Annals of
Statistics, 42(1):324–351, 2014.
[4] Laszlo Gyrofi, Michael Kohler, Adam Krzyzak, and Harro Walk. A Distribution-Free Theory of
Nonparametric Regression. 2002.
[5] David Heckerman, David Maxwell Chickering, David Meek, Robert Rounthwaite, and Carl
Kadie. Dependency networks for inference, collaborative filtering, and data visualization. Journalof Machine Learning Research, 1:49–75, 2000.
[6] Tao Hong, Pierre Pinson, and Shu Fan. Global energy forecasting competition 2012. International
Journal of Forecasting, 30:357–363, 2014.
[7] Martin Wainwright. Sharp thresholds for high-dimensional and noisy sparsity recovery usingconstrained quadratic programming (Lasso). IEEE Transactions on Information Theory, 55(5):2183–2202, 2009.
[8] Matt Wytock and Zico Kolter. Sparse Gaussian conditional random fields: Algorithms, theory,
and application to energy forecasting. In Proceedings of the 30th International Conference on
Machine Learning, pages 1265–1273, 2013.

16