Model selection and local geometry.
Robin J. Evans
University of Oxford
evans@stats.ox.ac.uk

arXiv:1801.08364v2 [] 31 Jan 2018

February 2, 2018
Abstract
We consider problems in model selection caused by the geometry of models close
to their points of intersection. In some cases, including common classes of causal
or graphical models as well as time series models, distinct models may nevertheless
have identical tangent spaces. This has two immediate consequences: first, in order to
obtain constant power to reject one model in favour of another we need local alternative
hypotheses that decrease to the null at a slower rate than the usual parametric n−1/2
(typically we will require n−1/4 or slower); in other words, to distinguish between the
models we need large effect sizes or very large sample sizes. Second, we show that
under even weaker conditions on their tangent cones, models in these classes cannot be
made simultaneously convex by a reparameterization.
This shows that Bayesian network models, amongst others, cannot be learned directly with a convex method similar to the graphical lasso. However, we are able to
use our results to suggest methods for model selection that learn the tangent space
directly, rather than the model itself. In particular, we give a generic algorithm for
learning discrete ancestral graph models, which includes Bayesian network models as
a special case.

1

Introduction

Consider a class of probabilistic models Mi indexed by some set i ∈ I, and suppose that we
have data from some distribution P ; model selection is the task of deducing, from the data,
which Mi contains P . Typically there will be multiple such models, in which case one may
appeal to parsimony or (if the model class is closed under intersection) select the smallest
such model by inclusion.
There have been dramatic advancements in certain kinds of statistical model selection,
including methods for working with large datasets and very high-dimensional problems (see,
for example, Bühlmann and van de Geer, 2011). However, model selection in some settings
is more difficult; for example, selecting an optimal Bayesian network for discrete data is
known to be an NP-complete problem (Chickering, 1996). In this paper we consider why
some model classes are so much harder to learn with than others. Taking a geometric
approach, we find that some classes contain models which are distinct but—in a sense that
will be made precise—are locally very similar to one another. The task of distinguishing
between them using data is therefore fundamentally more difficult, both statistically and
computationally.
Example 1.1. To illustrate the main idea in simple terms, consider a model space smoothly
described by a two dimensional parameter θ = (θ1 , θ2 )T , and with four submodels of interest:
M∅ : θ1 = θ2 = 0

M1 : θ2 = 0

M2 : θ1 = 0

M12 : unrestricted.
1

O(θ1 )

O(θ12 )

θ1

θ1

(a)

(b)

Figure 1: Illustration of model selection close to points of intersection. In (a) the models
M1 and M2 (in blue and red respectively) have different tangent spaces, so the distance
between them increases linearly as one moves away from the intersection. On the right, the
models M1 (in blue) and M02 (in red) have the same tangent space at the intersection, so
they diverge only quadratically with distance from the intersection.

In a setting with independent data, we would expect to have statistical power sufficient
to distinguish between M12 and M2 (i.e. to determine whether or not θ1 = 0) provided
that the magnitude of θ1 is large compared to n−1/2 , where n is the number of samples
available. We might also expect to be able to distinguish between M1 and M2 at the same
asymptotic rate; this is the picture in Figure 1(a), in which the distance between the two
models is proportional to distance from their intersection M∅ = M1 ∩ M2 (the constant of
proportionality being determined by the angle between the two models).
Suppose now that we define a model M02 : ψ1 = 0, where ψ1 ≡ θ12 − θ2 , and have to select
between M∅ , M1 , M02 , M12 (note that we still have M∅ = M1 ∩ M02 ), as illustrated in
Figure 1(b). Superficially, the task of choosing between these four models seems no different
to our first scenario, but in fact the models M1 and M02 are locally linearly identical at the
point of intersection ψ1 = θ2 = 0: that is, the tangent spaces of the two models at this point
are the same, so up to a linear approximation they are indistinguishable.
Models that overlap linearly in the manner above lead to two major consequences relating
to statistical power and computational efficiency. First, as illustrated in Figure 1(b), if M1
is correct the distance1 between the true parameter value (θ1 , 0) and the closest point on
M02 grows quadratically rather than linearly in θ1 . Hence, while |θ1 | = Ω(n−1/2 ) is sufficient
to gain power against M∅ , one needs |θ1 | = Ω(n−1/4 ) to ensure power2 against M02 . This is
potentially a very stringent condition indeed: if the effect size is halved, then we will need
16 times the sample size to maintain power against the alternative model.
Second, if two models have the same tangent space then we cannot choose a parameterization under which both models are convex sets. Note that in Figure 2(a) all four models are
convex, but in Figure 2(b) the model M02 is not. If we reparameterize to make M02 convex,
then M1 will not be3 . This prevents penalized methods such as the lasso being used in a
computationally efficient way.
Example 1.2 (Directed Gaussian Graphical Models). A common class of models in which
the phenomenon described above occurs is Gaussian Bayesian networks. Consider the two
graphs shown in Figure 2, each representing certain multivariate Gaussian distributions over
variables X, Y, Z with joint correlation matrix Π. The graph in Figure 2(a) corresponds to
1 Typically this would be approximated by the Mahalanobis distance using the Fisher information. For
regular statistical models, this is locally equivalent to the Hellinger distance or the square-root of the KLdivergence.
2 Recall that for positive functions f, g we have f (x) = Ω(g(x)) if and only if g(x) = O(f (x)).
3 By ‘reparameterize’ we mean under a twice differentiable bijection with an invertible Jacobian, so we
do not allow the map (θ1 , θ2 ) 7→ (ψ1 , θ2 ), for example.

2

X

Y

X

Y

Z

Z

(a)

(b)

Figure 2: Two Bayesian networks in which, for Gaussian random variables, the tangent
spaces of the models are identical at some points of intersection.

the marginal independence model X ⊥
⊥ Y , so that there is a zero in the corresponding
entry in Π: ρxy = 0. Figure 2(b), on the other hand, corresponds to the conditional
independence model X ⊥
⊥ Y | Z; that is, to a zero in the X, Y entry of Π−1 , or equivalently
to ρxy − ρxz ρzy = 0.
The two models intersect along two further submodels: for example, if ρxz = 0 (so
that X ⊥
⊥ Z), then ρxy = 0 if and only if ρxy − ρxz ρzy = 0. The same thing happens if
ρyz = 0 (i.e. Y ⊥
⊥ Z). When we are at the intersection between all these submodels—so
ρxy = ρxz = ρyz = 0 and all variables are jointly independent—we find that the tangent
spaces of the two original models are the same, giving rise to the phenomenon described
above. Indeed, we will see that this arises whenever two models intersect along two or more
such—suitably distinct—further submodels (Theorem 3.2).
In this paper we develop the notion of using local geometry as a heuristic for how closely
related two models are, and how rich a class of models is. In several classes of models for
which model selection is known to be difficult, we find that they contain distinct models
that are linearly equivalent at certain points in the parameter space.

1.1

Background and Prior Work

Crudely speaking, there are two flavours of statistical model, and consequently two main
reasons for wishing to select one. The first, called substantive or explanatory, emphasizes
the use of models to explain underlying phenomena, and such models are sometimes viewed
as approximations to an unknown scientific ‘ground truth’ (Cox, 1990; Skrondal and RabeHesketh, 2004). The second kind, referred to as empirical or predictive, is mainly concerned
with predicting outcomes from future observations, generally assuming that such observations will arise from the same population as previous data (Breiman, 2001). A discussion of
these two camps, together with some finer distinctions can be found in Cox (1990).
Our focus will primarily be on substantive models, in which case different models may
lead to rather different practical conclusions, even if the probability distributions associated
with them are ‘close’ in the sense observed above. The case of causal models such as the
graphs in Figure 2 is particularly stark: the reversal of an arrow will significantly affect
our understanding of how a system will behave under an intervention. For interpolative
prediction—that is, with new data from the same population as the data used to learn the
model—such concerns are generally lessened: if two probability distributions are similar
then they should give similar conclusions. However, the computational concerns we raise
will affect model selection performed for whatever reason.
For Bayesian network (BN) models specifically, there has been a great deal of work
dealing with the problem of accurate learning from data. Chickering (1996) showed that the
problem of finding the BN which maximizes a penalized likelihood criterion is NP-complete
in the case of discrete data with several common penalties. Uhler et al. (2013) give geometric
proofs that directed Gaussian graphical models are, in a global sense, very hard to learn
using sequential independence tests; this is because the volume of ‘unfaithful’ points that
will mislead at least one hypothesis test for a given sample size is very large, and in p > n
settings will overwhelm the model. Our approach is considerably simpler and is applicable to
3

arbitrary model classes and model selection procedures, but cannot make global statements
about the model.
Shojaie and Michailidis (2010) provide a penalized method for learning sparse highdimensional graphs, but they assume a known topological ordering for the variables in the
graph: that is, the direction of each possible edge is known. Ni et al. (2015) develop a
Bayesian approach that is similar in spirit, and apply it to gene regulatory networks. Fu
and Zhou (2013), Gu et al. (2014) and Aragam and Zhou (2015) all use penalization to
learn BNs without a pre-specified topological order, in the former paper even allowing for
interventional data; however in each case the resulting optimization problem is non-convex.
Other approaches based on assumptions such as non-Gaussianity or non-linearity are also
available (Shimizu et al., 2006; Bühlmann et al., 2014).
The remainder of the paper is organized as follows. In Section 2 we formalize the intuition
given above by carefully defining local similarity between models, and then proving results
relating to local asymptotic power and convex parameterizations. In Section 3 we give
sufficient conditions for this situation to occur, including the intersection of two models
along multiple distinct submodels. In Section 4 we apply these results to show that the lasso
cannot be used directly to learn Bayesian networks. Section 5 provides further examples of
how the result can be applied, while Section 6 considers related phenomena such as double
robustness. Section 7 suggests methods to exploit model classes in which the tangent spaces
are distinct but in which we cannot make the models simultaneously convex, and Section 8
contains a discussion.

2

Models

Consider a class of finite-dimensional probability distributions {Pθ : θ ∈ Θ}, where Θ is an
open subset of Rk . We assume throughout that the parameter θ describes a smooth (twice
differentiable with full rank Jacobian) bijective map between the set of distributions and Θ.
Consequently, we will refer interchangeably to a subset of parameters and the corresponding
set of probability distributions as a model.
Suppose we have models corresponding to subsets Mi ⊆ Θ, i = 1, 2, . . .. Some of our
results are applicable specifically when these subsets are algebraic varieties—that is, those
defined by the zeroes of a set of polynomials. This includes a wide range of models of
interest, see Drton and Sullivant (2007) for details. Our formal discussion of the similarity
of these models is based partly on their tangent spaces at points of intersection.
Definition 2.1. The tangent cone Cθ (M) of a model M at θ ∈ M is defined as the set of
limits of sequences αn (θn − θ), such that αn > 0, θn ∈ M and θn → θ.
The tangent space of M at θ is the vector space Tθ (M) spanned by elements of the
tangent cone.
Clearly the tangent space contains the tangent cone; the model is said to be regular at
θ if the two are equal; in particular this means that M looks like a differentiable manifold
(embedded in Θ) at θ, and regular parametric asymptotics apply: in other words, the
maximum likelihood estimator is asymptotically normal with covariance given by the inverse
Fisher information (van der Vaart, 1998). Most of the models we will initially consider are
regular everywhere, though their intersections may not be.
Our next definition considers the classification of models based on their local similarity.
We work with a local version of the Hausdorff distance between sets; this is the furthest
distance from any point on one set to the nearest point on the other. Denote this by


D(A, B) ≡ max sup inf ka − bk, sup inf ka − bk .
b∈B a∈A

a∈A b∈B

4

Definition 2.2. We say that M1 , M2 are c-equivalent at θ ∈ M1 ∩ M2 if the Hausdorff
distance between the sets in a ball of radius ε is o(εc ). Formally:
lim ε−c D(M1 ∩ Nε (θ), M2 ∩ Nε (θ)) = 0,
ε↓0

where Nε (θ) is an ε-ball around θ. In other words, within an ε-ball of θ, the maximum
distance between the two models is o(εc ).
If the limit above is finite but not necessarily zero—i.e. the distance is O(εc )—we will
say that M1 and M2 are c-near-equivalent at θ.
The definition of c-equivalence is given by Ferrarotti et al. (2002), who also derive some
of its elementary properties. Given that θ ∈ M1 ∩ M2 then 1-near-equivalence is trivial,
while 1-equivalence means that to a linear approximation around θ the models are the same
(formally they have the same tangent cone at θ); this is illustrated by the two surfaces in
Figure 3(a). Similar considerations can be applied to higher orders: 2-equivalence means that
the quadratic surface best approximating one model is the same as that approximating the
other. For smooth models, (k−1)-equivalence for k ∈ N typically implies k-near-equivalence.
The idea of 1-equivalence is similar to what Le Cam (1960) (see also the overview by
van der Vaart, 2002) calls asymptotically differentially equivalent models, but that concept
applies at all parameter values, not just a specific one.
Proposition 2.3 (Ferrarotti et al. (2002), Proposition 1.3). Two models M1 , M2 are 1equivalent at θ ∈ M1 ∩ M2 if and only if Cθ (M1 ) = Cθ (M2 ).
Example 2.4. Consider the two graphical models in Figure 2, defined respectively by the
independence constraints X ⊥
⊥ Y and X ⊥
⊥ Y | Z. If we take the set of trivariate Gaussian
distributions satisfying these two restrictions then, as discussed in Section 1, the two models
are 2-near-equivalent at all diagonal covariance matrices.
For finite discrete variables, however, these two models are not 1-equivalent. Suppose
X, Y, Z take nX , nY , nZ levels respectively. The models are defined by the equations
X⊥
⊥Y |Z

:

p(x, y, z) · p(z) − p(x, z) · p(y, z) = 0

X⊥
⊥Y

:

p(x, y) − p(x) · p(y) = 0

∀x, y, z
∀x, y,

where, for example, p(x, z) = P (X = x, Z = z). In particular, the marginal independence
model is subject to (nX − 1)(nY − 1) restrictions, but the conditional independence model
is subject to that many restrictions for each level of Z. One can use this to show that the
dimension of the conditional independence model is smaller than the marginal independence
model, and so they cannot both be approximated by the same linear space.
In spite of the differing dimensions in the example above, we will show that the issue of
the models not being simultaneously convex will still arise. This motivates a slightly weaker
definition for models meeting in an intuitively ‘irregular’ manner, which we term overlap.
Definition 2.5. Suppose we have two models M1 , M2 with a common point θ ∈ M1 ∩M2 .
We say that M1 and M2 overlap if there exist points h ∈ Cθ (M1 ) ∩ Cθ (M2 ) but h 6∈
Cθ (M1 ∩ M2 ).
In other words, there are tangent vectors that can be obtained in either model, but not
in their intersection. Note that necessarily we have Cθ (M1 ∩ M2 ) ⊆ Cθ (M1 ) ∩ Cθ (M2 ),
so the definition asks that this is a strict inclusion. For distinct, regular algebraic models,
overlap is implied by 1-equivalence.
Lemma 2.6. Let M1 and M2 be algebraic varieties that are 1-equivalent and regular at θ,
but not equal in a neighbourhood of θ. Then M1 and M2 overlap at θ.

5

(a)

(b)

Figure 3: Illustration of two models that (a) are 1-equivalent; (b) overlap. In (a) the two
surfaces have the same tangent space, whereas in (b) they do not.
Proof. Suppose that, at θ, the models do not overlap and are 1-equivalent. We will show
that they are equal in a neighbourhood of θ.
Since M1 is regular at θ we can assume that it is an irreducible variety, else replace it with
its unique irreducible component containing θ (similarly for M2 ). The condition that the
models are 1-equivalent implies Cθ (M1 ) = Cθ (M2 ), and no overlap means Cθ (M1 ∩ M2 ) =
Cθ (M1 ) ∩ Cθ (M2 ) = Cθ (M1 ).
On the other hand, M1 ∩ M2 is a subvariety of M1 , so it is either equal to M1 or
has strictly lower dimension. If the latter, then Cθ (M1 ∩ M2 ) would be of this same lower
dimension (Cox et al., 2008, Theorem 9.7.8), and hence is a strict subset. Since we have
already seen that M1 ∩ M2 has the same tangent cone as M1 , it must be that the former
possibility holds; that is, M1 ∩ M2 = M1 , in a neighbourhood of θ. Clearly the same is
true for M2 , which proves the result.
In fact, it seems likely that a substantially weaker condition than the models being
algebraic would suffice, but it is not clear to us exactly what that condition should be.
A canonical example of sets that overlap but are not 1-equivalent is given by the subsets
of R3 shown in Figure 3(b); M1 = {(x, y, z) : y = z = 0} (in blue) and M2 = {(x, y, z) : z =
−x2 } (in red) have M1 ∩ M2 = {0}. The tangent cone of M2 is the plane z = 0, while M1
is its own tangent cone and therefore C0 (M1 ) ⊆ C0 (M2 ) and C0 (M1 ) ∩ C0 (M2 ) = M1 .
However, C0 (M1 ∩ M2 ) = {0}, so the inclusion is strict. In words, we can approach the
origin along a line tangent to the x-axis in either model, but not in the intersection. The
blue model has smaller dimension than the red so they clearly not 1-equivalent.
Typically it is hard to show that two models are not 1-equivalent (or do not overlap)
anywhere in the parameter space: indeed the difficulty of verifying such global conditions
is one of our motivations for considering these local criteria instead. However, for algebraic
models we can often show that such points are at most a set of zero measure within a large
class of ‘interesting’ submodels.

2.1

Statistical Power

We will assume that all the models we consider satisfy standard parametric regularity conditions, in particular differentiability in quadratic mean (DQM), which yields the familiar
asymptotic expansion of the log-likelihood `(θ) (see, e.g. van der Vaart, 1998):
hT ˙
1
`(θ + hn ) − `(θ) = √ `(θ)
− hT I(θ)h + op (1),
2
n

(1)

˙
where `(θ)
is the data dependent score, I(θ) the Fisher information for one observation,
and n1/2 hn → h. For the purpose of distinguishing between models from data we need the
6

difference between the log-likelihoods at points close to the MLE not to vanish as sample
size n → ∞. The expansion above shows that this requires h 6= 0 for the right hand side to
contain a stable term. Hence the distance between the two parameter values needs to shrink
no faster than n−1/2 , the standard parametric rate of statistical convergence.
We consider settings in which hn , h̃n → 0 with n1/2 (hn − h̃n ) → k, so that we may
compare alternatives in two different models. We also impose the mildly stronger condition
that the model is uniformly differentiable in quadratic mean (UDQM) in some neighbourhood of θ (see Appendix A.2 for details). This will hold, for example, on the interior of the
parameter space of any exponential family model.
Theorem 2.7. Let M be a model that is UDQM at some θ, and M1 , M2 ⊂ M be submodels
with θ ∈ M1 ∩ M2 . Let hn , h̃n → 0 be sequences such that θ + hn ∈ M1 and θ + h̃n ∈ M2
and k = limn n1/2 (hn − h̃n ). Then
1
kT ˙
+ h̃n ) − k T I(θ)k + op (1)
`(θ + hn ) − `(θ + h̃n ) = √ `(θ
2
n


1
d
T
T
−→ N − k I(θ)k, k I(θ)k .
2
The proof is given in Appendix A.2. The theorem shows that it is not possible to
distinguish between the two models as n grows if k = 0, that is, if the difference in the
alternatives shrinks at a rate faster than n−1/2 . The geometry of the model determines
the relationship between the rate at which hn , h̃n shrink and the size of h̃n − h; in general
kh̃n − hk may be of smaller order than both khn k and kh̃n k.
This leads to our first main result about c-equivalence, which is a direct application of
its definition.
1

− 2c
) with
Theorem 2.8. Suppose M1 , M2 are c-equivalent at θ, and
√ consider hn = O(n
θ + hn ∈ M1 . There exists h̃n with θ + h̃n ∈ M2 such that n(hn − h̃n ) → 0.

Combining the previous two theorems gives the corollary of interest.
Corollary 2.9. If two models are c-equivalent (repectively c-near equivalent) at θ then
1
they cannot be asymptotically distinguished under local alternatives (to θ) of order O(n− 2c )
1
(respectively o(n− 2c )).
Most models are 1-near-equivalent, and therefore we recover the usual parametric rate:
we have power only if hn shrinks at n−1/2 or slower. For 1-equivalent models this is not
enough: a rate of n−1/2 will be too fast for us to tell whether our parameters are in M1
or M2 . In practice, regular models that are 1-equivalent are also 2-near-equivalent, so any
rate quicker than n−1/4 is also too fast.

2.2

Convexity

A second consequence of having 1-equivalent regular models is that it is not possible to have
a parameterization with respect to which both models are convex; in fact, this is true even
under the weaker assumption of overlap. Some automatic model selection methods such as
the lasso (Tibshirani, 1996) rely on a convex parameterization, usually constructed by ensuring that interesting submodels correspond to coordinate zeroes (e.g. θ1 = 0). Such selection
methods cannot be directly applied to model classes that contain overlapping models.
Theorem 2.10. Suppose M1 and M2 are distinct models that overlap at a point θ on their
relative interiors. Then there is no parameterization with respect to which both these models
are convex in a neighbourhood of θ.

7

⇒

(a)

⇒

(b)

Figure 4: Cartoon illustrating why overlapping models cannot be made locally convex. In
(a) the red and blue models intersect at a non-zero angle and so the parameterization can
be smoothly transformed to make the models locally linear. In (b) there is no angle between
the models, and this fact is invariant in any smooth reparameterization.
Proof. This is just a statement about sets under differentiable maps with differentiable
inverses. We will prove the contrapositive: if two sets C, D with x ∈ C ∩ D are convex (in a
neighbourhood of x) then they do not overlap at x. Overlap is a property of tangent cones,
and these are structurally preserved under maps that are differentiable with an invertible
Jacobian at x; this will prove the result.
The affine hull of any convex set C is the unique minimal affine space A containing C.
Furthermore, by definition of the relative interior, for every x ∈ C int there is a neighbourhood
of x in which C and A are identical, i.e. C coincides with the affine space A at x.
Now, suppose two sets C and D are both convex and x ∈ C int ∩ Dint . Then they coincide
with affine spaces say A, B and hence are their own tangent spaces at x. This implies that
C ∩ D coincides with the affine space A ∩ B, and hence the tangent cone of C ∩ D is just
A ∩ B. Hence C and D do not overlap.
Intuitively, one can reparameterize a regular model such that (at least locally to some
point θ) the model is convex. However, this cannot be done simultaneously for two overlapping models. This is illustrated in Figure 4. The result fails if θ is not required to
be a point on the relative interior, a counterexample being models M1 = {θ2 ≥ θ12 } and
M2 = {θ2 ≤ −θ12 } with θ = (0, 0).
Remark 2.11. Many convex methods, including the lasso, proceed by optimizing a convex
function over a parameter space containing all submodels of interest. In the case of the lasso,
this function has ‘cusps’ on submodels of interest that lead to a non-zero probability that
the optimum lies exactly on the model. Theorem 2.10 shows that such a procedure cannot
be convex if the class contains overlapping models, since the submodels themselves cannot
be made convex. Any submodel that is not convex cannot represent a cusp in a convex
function, and therefore we cannot obtain a non-zero probability of selecting such a model.
The nature of this impossibility result should perhaps not be surprising, since all versions
of the lasso in the context of ordinary linear models place a restriction on the collinearity
of the different parameters in the form of restricted isometry properties or similar (for an
overview see, for example, Bühlmann and van de Geer, 2011). In the case of 1-equivalent
models, as we move closer to the point of intersection the angle between the two models
shrinks to zero, so no such property could possibly hold. For models that overlap the problem
is similar, but may only apply as we approach from certain directions. The result implies,
in particular, that an algorithm such as the lasso cannot be used directly to learn Bayesian
networks, whether Gaussian or discrete. We expand upon this in Section 4.
If a class of models is non-convex when parameterized in a canonical way, it may be
possible to reparameterize so that they are all convex, but only if no two models overlap.
8

For example, suppose we are interested in Gaussian models of marginal independence, that
is, defined by the pattern of zeroes in the off-diagonal elements of the covariance matrix.
The log-likelihood of a multivariate Gaussian with zero mean and covariance matrix Σ is
`(Σ; S) =

n
log det Σ−1 − tr(SΣ−1 ) ,
2

where S is the sample covariance matrix and n is the sample size. This is a simple function
of the inverse covariance matrix Σ−1 —the canonical parameter for this exponential family—
but a complicated function of the mean parameter Σ. Indeed, while the log-likelihood is a
convex function of Σ−1 , it is typically not as a function of Σ, and placing linear restrictions
on this parameter leads to problems with multiple local maxima (see Zwiernik et al., 2017,
for an overview of related problems).
However, if we are prepared to accept some loss of efficiency, there is nothing to stop us
estimating Σ via a moment matching approach, such as by solving the convex program:




X
Σ̂ = arg minΣ0 kΣ − Sk2 + ν
|σij | ;


i<j

here Σ  0 denotes that Σ belongs to the convex set of positive definite symmetric matrices.
If the penalty ν is chosen to grow at an appropriate rate (nδ for some 12 < δ < 1) then
under some conditions on the Fisher information for Σ, this will still be consistent for model
selection (see, for example, Rocha et al., 2009, Theorem 5).

3

Submodels, Equivalence and Overlap

In this section we consider sufficient conditions for models to be c-equivalent or to overlap.
We will see that in many cases it is a consequence of having models whose intersections are
themselves expressible as a union of two or more distinct models.

3.1

c-equivalence

For the example in Figure 2, the 1-equivalence of the two models M1 = {Π : X ⊥
⊥ Y } and
M2 = {Π : X ⊥
⊥ Y | Z} was closely related to the fact that if either X ⊥
⊥ Z or Y ⊥
⊥Z
holds, then M1 and M2 intersect. This means that their intersection is a union model ;
that is, it can be non-trivially expressed as the union of two or more models. It follows
that at any points where both the submodels MX⊥
⊥Z and MY ⊥
⊥Z hold, the two original
models are 1-equivalent. This is because any direction in M1 (or M2 ) away from the point
of intersection M1 ∩ M2 can be written as a linear combination of (limits of) vectors that
lie in one of MX⊥
⊥Z or MY ⊥
⊥Z ; the partial derivatives of the distance between M1 and
M2 in these directions are zero, and so any directional derivative of this quantity is zero.
This is illustrated in Figure 5(a), which shows two surfaces intersecting along two lines: in
directions that are diagonal to these lines, the two surfaces separate at a rate that is at most
quadratic.
It is not necessary for models to intersect in this manner in order to be 1-equivalent:
for example, the surfaces z = 0 and z = x2 + y 2 intersect only at the point (0, 0, 0) (see
Figure 5(b)). However, many models that are 1-equivalent do intersect along at least two
submodels, including most of the substantive examples that we are aware of; the time series
models in Section 5.2 are an exception. If k > 2 submodels are involved in the intersection,
then we will see that the original models are k-near-equivalent and asymptotic rates for local
alternatives will be even slower than n−1/4 .
The next Proposition provides a basis for proving that result, as well as giving some
intuition in terms of Taylor expansions.

9

(b)
(a)

Figure 5: Illustration of two surfaces that are 1-equivalent such that: (a) they intersect
along two submodels; (b) they do not.
Proposition 3.1. Let f : X ⊂ Rk → R be a C m function such that f (0) = 0, and let hji :
X → R, i = 1, . . . , m, j = 1, . . . , ci be a collection of C m functions such that hji (0) = 0, ∀i, j
and the Jacobian J(h11 , . . . , hcmm ) has full rank at x = 0. Suppose further that f (x) = 0
whenever h1i (x) = · · · = hci i (x) = 0 for any i ∈ {1, . . . , m} and x ∈ X. Then f (x) =
O(kxkm ).
Proof. First, reparameterize the space X so that hji = 0 corresponds to a coordinate zero.
Since all the models are C m and using the Jacobian condition, this can be done with a C m
function.
Now, since f is an m-times differentiable function we can write out its Taylor expansion
at 0 up to order m and have a remainder of order o(kxkm ). In order for f to be non-zero,
we need to move off the submodels defined by h1i = · · · = hci i = 0 for every i = 1, . . . , m.
The first non-zero partial derivatives will then be of the form
∂mf
,
∂z1 · · · ∂zm
where there is one zi normal to each submodel.
This phenomenon is described in the algebraic geometry literature as a zero of multiplicity
m. It follows from this result that if we have two models that intersect along two distinct
submodels then they will be 2-near-equivalent; if they intersect along three or more distinct
submodels, then they share not only a tangent space but also higher-order features such as
curvature.
To formalize this, we use the next result. Define the normal space of M at θ, denoted
Tθ (M)⊥ , to be the orthogonal complement to the tangent space (within some suitably
defined ambient space).
Theorem 3.2. Suppose M1 , M2 and N1 , . . . , Nm are C m manifolds all containing a point
θ, and such that Ni ∩ M1 = Ni ∩ M2 for each i. Suppose further that the normal spaces of
N1 , . . . , Nm at θ are all linearly independent. Then M1 and M2 are m-near-equivalent at
θ.
In other words, there are m distinct submodels on which M1 and M2 intersect; in order
to distinguish between M1 and M2 , we need to ‘move away’ from all the submodels Ni .
The condition on the normal spaces is to ensure that we cannot move directly away from
several submodels at once.
Proof. See Appendix.

10

1

2

3

1

2

3

4

4

(a)

1

(b)

2

...

k−1

k

k+1
Figure 6: (a) and (b) are two graphs that differ only by the arrow heads present at 3. (c) A
discriminating path of length k.

Perhaps an easier way to think about the normal spaces is in terms of the submodels
being defined by ‘independent constraints’; in particular, by constraints defined on different
parts of the model. If N1 is defined by the set of points that are zeros of the functions
f1 , . . . , fk , then Tθ (N1 )⊥ contains the space spanned by the Jacobian J(f1 , . . . , fk ); if this
is of full rank, then the two spaces are equal. If N2 is similarly defined by g1 , . . . , gl then the
condition is equivalent to saying that the Jacobian of all k + l functions has full rank k + l.
Example 3.3 (Discriminating Paths). The graphs in Figures 6(a) and (b) are examples
of ancestral graphs, which will be introduced more fully in Section 4. Both these graphs
are associated with probabilistic models in which X1 ⊥
⊥ X3 , but in the case of (a) we also
have X1 ⊥
⊥ X4 | X2 , X3 , whereas (b) implies X1 ⊥
⊥ X4 | X2 . The difference between
the conditioning sets in the latter two independences makes the models distinguishable.
However, if any of the three edges between the pairs (1,2), (2,3) and (3,4) are removed,
then the models become Markov equivalent—that is, they represent the same conditional
independences and therefore are identical. In the case of multivariate Gaussian distributions
this corresponds to any of the (partial) correlations ρ12 , ρ23 or ρ34·12 being zero. Applying
Theorem 3.2 shows that these models are 3-near-equivalent at points where all variables are
independent.
This example can be expanded to arbitrarily long discriminating paths of the kind shown
in Figure 6(c): these differ only by the edges incident to the vertex k, and this makes them
distinguishable. However, if any of the edges (i, i + 1) for i = 1, . . . , k are missing, the
submodels are Markov equivalent, and hence these Gaussian graphical models are k-nearequivalent. This means that the ‘discrimination’ between different models that is possible
may be quite limited in practice, absent extraordinarily large sample sizes. We provide a
simulation study to illustrate this in Section 4.4.
This example has serious ramifications for the FCI algorithm, which uses discriminating
paths to orient edges (Zhang, 2008). It suggests only very strong dependence will allow an
unambiguous conclusion to be reached for moderate to long paths.

3.2

Identifying Overlap

Overlap between two regular models occurs when their intersection is not itself a regular
model, but rather a union of such models. Let M be an algebraic model that can be written
as M = V1 ∪ V2 for incomparable varieties V1 and V2 ; in this case we say M is reducible,
and otherwise irreducible.

11

Theorem 3.4. Let M1 and M2 be two regular algebraic models, and suppose that M1 ∩M2
is reducible into two (or more) further models that have incomparable tangent cones at
θ ∈ M1 ∩ M2 . Then the models M1 and M2 overlap at θ.
This condition is similar in spirit to Theorem 3.2, but note that here we do not require the
normal spaces of the two submodels to be linearly independent at the point of intersection,
only that their tangent spaces are incomparable.
Proof. Let M1 ∩ M2 = V1 ∪ V2 be the reduction into subvarieties. Taking θ ∈ V1 ∩ V2 , note
that Cθ (M1 ∩ M2 ) = Cθ (V1 ∪ V2 ) = Cθ (V1 ) ∪ Cθ (V2 ), the second equality following from the
definition of a tangent cone. Now, since V1 , V2 ⊆ M1 and M1 is regular at θ, this implies
that any vector in Cθ (V1 ) + Cθ (V2 ) is contained in Cθ (M1 ); similarly for M2 . Therefore
the condition for not overlapping,
Cθ (M1 ∩ M2 ) = Cθ (M1 ) ∩ Cθ (M2 ),
holds only if Cθ (V1 ) + Cθ (V2 ) ⊆ Cθ (V1 ) ∪ Cθ (V2 ), which occurs only if one of Cθ (V1 ) or
Cθ (V2 ) is a subspace of the other; this was ruled out by hypothesis.
Example 3.5. As already noted in Examples 1.2 and 2.4, the Gaussian graphical models
defined respectively by the independences M1 : X ⊥
⊥ Y and M2 : X ⊥
⊥ Y | Z are 1equivalent at diagonal covariance matrices, but the corresponding discrete models are not.
This is because—taking X, Y, Z to be binary—the three-way interaction parameter
X
λXY Z ≡
(−1)|x+y+z| log P (X = x, Y = y, Z = z)
x,y,z∈{0,1}

is zero in the conditional independence model, but unrestricted in the marginal independence
model. However, the intersection of M1 and M2 for binary X, Y, Z is the set of distributions
such that either X ⊥
⊥ Y, Z or Y ⊥
⊥ X, Z, and these correspond respectively to the submodels
λXY = λXZ = λXY Z = 0

or

λXY = λY Z = λXY Z = 0,

defined by zeros of log-linear parameters (see Appendix A.3 for full definitions). These
models satisfy the conditions of Theorem 3.4, and therefore M1 and M2 overlap at points
of total independence ⊥
⊥{X, Y, Z}.

4

Directed and Ancestral Graph Models

In this section we focus on two classes of graphical models: Bayesian network models, and
the more general ancestral graph models. A more detailed explanation of the relevant theory
can be found in Spirtes et al. (2000) and Richardson and Spirtes (2002).

4.1

Ancestral Graphs

A maximal ancestral graph (MAG) is a simple, mixed graph with three kinds of edge,
undirected (−), directed (→) and bidirected (↔). Special cases of ancestral graphs include
directed acyclic graphs, undirected graphs and bidirected graphs, but not chain graphs.
There are some technical restrictions on the structure of the graph which we omit here
for brevity: the key detail is that—under the usual Markov property—the model implies
a conditional independence constraint between each pair of vertices if (and only if) they
are not joined by any sort of edge in the graph (Richardson and Spirtes, 2002). The set
that needs to be conditioned upon to obtain the independence depends on the presence of
colliders in the graph. A collider is a pair of edges that meet with two arrowheads at a
vertex k: for example, i → k ← j or i → k ↔ j. We say the collider is unshielded if i and j
are not joined by an edge. Any other configuration is called a noncollider.
12

The special case of an ancestral graph model in which all edges are directed yields a
Bayesian network (BN) model, widely used in causal inference and in machine learning
(Bishop, 2007; Pearl, 2009). The additional undirected and bidirected edges allow MAGs
to represent the set of conditional independence models generated by marginalizing and
conditioning a BN model. Ancestral graphs are therefore useful in causal modelling, since
they represent the conditional independence model implied by a causal structure with hidden
and selection variables.
Example 4.1. Consider the maximal ancestral graphs in Figure 6(a) and (b). The graph
in (a) is fully directed and represents the model defined by the conditional independences:
X1 ⊥
⊥ X3 ,

X1 ⊥
⊥ X4 | X2 , X3 .

The graph in (b), on the other hand, represents
X1 ⊥
⊥ X3 ,

X1 ⊥
⊥ X4 | X2 .

The difference in the conditioning sets above is due to the fact that 2 ← 3 → 4 is a noncollider
in the first graph, but a collider in the second: 2 ↔ 3 ↔ 4.
An independence model, I, is a collection of (conditional) independence statements of
the form Xi ⊥
⊥ Xj | XC , for i 6= j and possibly empty C. We will say I is simple if it
contains at most one independence statement for each unordered pair {i, j}. If there is no
independence statement between Xi , Xj in a simple independence model, we say i and j are
adjacent.
Theorem 4.2. Let I1 , I2 be simple independence models on the space of p × p Gaussian
covariances matrices. Then the two models are 2-near-equivalent if they have the same
adjacencies.
Further, if the two models have different adjacencies then they are 1-equivalent on at
most a null set within any parametric independence model.
Proof. Let E denote the set of adjacencies in a simple independence model I. Parameterizing
using the set of correlation matrices, we will show that the tangent space of I at the identity
matrix I is
M
TI (I) =
Dij ,
(2)
i<j
{i,j}∈E

where the matrices Dij have zeroes everywhere except in the (i, j) and (j, i)th entries, which
are 1.
If {i, j} ∈ E it is easy to see that I + λDij is in the model for all λ ∈ (−1, 1), since
all conditional independences except those between Xi and Xj hold; hence Dij ∈ TI (I).
Conversely, if i and j are not adjacent, then some independence restriction Xi ⊥
⊥ Xj | XC
holds, so
f (Π) = Πij − ΠiC (ΠCC )−1 ΠCj = 0.
The derivative of f at Π = I is just Dij , so it follows that Dij 6∈ TI (I). Hence the tangent
space at I is in the form (2). By Proposition 2.3, the models are 1-equivalent, and since
these constraints are linearly independent at Π = I, they are regular and therefore also
2-near-equivalent.
Conversely, suppose that there is some pair i, j subject to the restriction Xi ⊥
⊥ Xj | XC
in I1 but not to any such restriction in I2 . By the above analysis, these models have
distinct tangent spaces at Π = I. Since these are models defined by polynomials in Π, the
set of points on which the tangent spaces are identical (say W ) is an algebraic variety; its
intersection with any irreducible variety V is therefore either equal to V or of strictly smaller
13

dimension than V (indeed this follows from the usual definition of dimension in such sets;
see Bochnak et al. (2013, Section 2.8)). However, if the identity matrix is contained in V
then clearly W ∩ V ⊂ V , since we have established that the tangent spaces do not intersect
at the identity. Hence W ∩ V has smaller dimension, and is a null subset of V .
Corollary 4.3. Two Gaussian maximal ancestral graph models are 2-near-equivalent if
they have the same adjacencies, and are 1-equivalent almost nowhere on any submodel of
independence otherwise.
Proof. Follows from the pairwise Markov property of Richardson and Spirtes (2002).
Corollary 4.4. Two Gaussian chain graph models, under any of the interpretations given
in Drton (2009), are 2-near-equivalent if they have the same adjacencies.
Proof. See Remark 5 of Drton (2009), which makes clear that pairwise independence are
sufficient to define Gaussian chain graph models.
We conjecture that, in fact, two ancestral graph models of the form given in Theorem
4.2 will overlap nowhere in the set of positive definite correlation matrices if they do not
share the same adjacencies (rather than almost nowhere). It is not hard to see that this
holds for models of different dimension, since these models are regular and therefore will
share this dimension everywhere. To prove it in general seems challenging; the result above
is sufficient for most practical purposes.

4.2

Discrete Data

The picture is slightly rosier if we consider discrete data instead.
A well known result of Chickering (1996) shows that finding an optimal Bayesian Network
for discrete data is an NP-hard problem; in other words, it is computationally difficult.
However, we find that from a statistical point of view, it is somewhat easier than in the
Gaussian case.
Theorem 4.5. No two distinct, discrete, maximal ancestral graph models are 1-equivalent
at an independence model.
Proof. We provide a sketch proof for binary variables; the general discrete case follows
similarly. Let G be a graph with vertices V . We can parameterize the binary probability
simplex using log-linear parameters λA for ∅ 6= A ⊆ V . We will consider the tangent space
of the model at the uniform distribution; that is, at λA = 0 for every A.
The conditional independence Xa ⊥
⊥ Xb | XC is equivalent to the λ0abD = 0 for each
0
D ⊆ C, where λA are the log-linear parameters for the marginal distribution over Xa , Xb , XC
(Rudas et al., 2010). However, within δ of the uniform distribution we have λ0A = λA +O(δ 2 )
(see Corollary A.5). Hence the constraint to the tangent cone imposed by λ0abD = 0 is the
same as that imposed by λabD = 0 for each D ⊆ C.
Now, two MAGs are Markov equivalent if and only if they have the same adjacencies,
unshielded colliders, and discriminating paths (e.g. Zhang, 2008, Proposition 2). If they
differ in adjacencies (say i, j), then a log-linear parameter λij will appear in the tangent
cone of one model but not the other. If they differ in an unshielded collider i ∗→ k ←∗ j
then in one model λijk = 0 but in the other this direction is not restricted. The same holds
for a discriminating path between i and j for a potential collider k.
Remark 4.6. Note that this proof also demonstrates that the skeleton of the graph is
determined solely by the two-way interaction parameters, and that the remainder of the
model can be deduced entirely from the three-way interaction parameters. This suggests
that it might be possible to develop a model selection procedure using only this information,
something we do in Section 7. It also illustrates that in cases with strong three-way interactions it should be easier to learn the correct model rather than just the correct skeleton. A
14

X

Y

X

Y

Z

W

Z

W

(a)

(b)

Figure 7: Two Bayesian networks which do not overlap but are not consistent with a single
topological order.

‘noisy-OR’ model would, for example, have the desired property. The well known ALARM
dataset (Beinlich et al., 1989) has strong interaction effects and is—at least in part for this
reason—considered to be relatively easy to learn.

4.3

Overlap

Although distinct, discrete MAG models (and therefore BN models) are never 1-equivalent,
they do still overlap, as our next result demonstrates.
Proposition 4.7. Let M(G1 ), M(G2 ) be two discrete Bayesian network models such that
i − k − j is an unshielded collider in G1 and an unshielded noncollider in G2 . Then the two
models overlap.
Proof. We claim that the space spanned by λik and λjk is contained in the tangent cones
of both models, but not of their intersection; the first claim follows from Theorem 4.5. For
the second, if Xk is binary then λij
ij = λij = 0 if and only if either λik = 0 or λjk = 0 (see
Drton et al., 2008, Example 3.1.7). Clearly then directions in which they are both non-zero
will not appear in the intersection model.
If Xk has more than two states then a similar conclusion can be reached.
Bayesian network models that are consistent with a single topological ordering of the
vertices do not overlap, because their intersection is always another BN model. We can
therefore work with a class defined by the subgraphs of a single complete BN in order to
avoid the problems associated with overlap.
This leads to the question of whether any other, perhaps larger, subclasses share this
property. The previous result shows that any such subclass would be restricted fairly
severely, since any two graphs must never disagree about a specific unshielded collider.
Note that the result does not imply that it is necessary for graphs to be consistent with a
single topological order in order for the corresponding models not to overlap; the graphs in
Figure 7 provide a counterexample to this.
The easiest way to ensure that a class of models does not overlap at the independence
model is to associate each potentially missing edge with a single constraint using a pairwise
Markov property, as with the set of BN models that are consistent with a given topological
order, or the set of undirected graph models.
Proposition 4.7 and Theorem 2.10 show that the impossibility result discussed in Remark
2.11 applies to Bayesian networks, and we will never be able to use a lasso-like method to
consistently select from this class under standard conditions.

4.4

Discriminating Paths

Example 3.3 introduced the notion of a discriminating path, which allows the identification
of colliders in an ancestral graph. Formally, define the ancestral graphs Gk and Gk0 as having
15

s
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5

ρs
0.4
0.283
0.2
0.141
0.1
0.071
0.05
0.035
0.025
0.018
0.012

k=2
n
32
128
512
2048
8 192
32 768
131 072
524 288
2 097 152
8 388 608
33 554 432

acc.
0.703
0.681
0.679
0.677
0.660
0.681
0.654
0.674
0.659
0.684
0.672

k=3
n
250
2 000
16 000
128 000
1 024 000
8 192 000
65 536 000
524 288 000
4.19 × 109
3.36 × 1010
2.68 × 1011

acc.
0.776
0.730
0.711
0.707
0.693
0.695
0.699
0.690
0.704
0.682
0.691

k=4
n
800
12 800
204 800
3 276 800
52 428 800
838 860 800
1.34 × 1010
2.15 × 1011
3.44 × 1012

acc.
0.763
0.669
0.639
0.651
0.652
0.645
0.628
0.650
0.652

k=5
n
5 000
160 000
5 120 000
163 840 000
5.24 × 109
1.68 × 1011
5.37 × 1012

Table 1: Table giving effect sizes ρs and sample sizes used in each of the graphs Gk for
k = 2, 3, 4, 5, and corresponding to the solid lines in Figure 8. For k = 4, 5 the sample sizes
grow extremely quickly, and some larger entries are excluded because they led to numerical
problems. ‘acc.’ gives the proportion of runs for which the correct graph was selected.
vertices 1, . . . , k+1, with a path 1 ↔ 2 ↔ · · · ↔ k and directed edges from each of 2, . . . , k−1
to k + 1. In addition, Gk has the edges k ↔ k + 1, while Gk0 has k → k + 1; the graphs are
shown in Figure 6(c), with only the final edge left ambiguous. The path from 1 to k + 1 is
known as a discriminating path, and its structure allows us to determine whether there is a
collider at k (as in Gk ) or not (as in Gk0 ).
As noted already, if any of the edges i ↔ i + 1 for i = 1, . . . , k − 1 are missing (i.e., if
Xi ⊥
⊥ Xi+1 ), then the two models coincide. In addition, if the final edge between k and
k + 1 is missing (so Xk+1 ⊥
⊥ Xk | X1 , . . . , Xk−1 ), the two models coincide. We consider the
Gaussian graphical models associated with these graphs and, by application of Theorem 3.2,
the two ancestral graph models for Gk and Gk0 are k-near-equivalent at points where these
additional independences hold.
We now perform a small simulation to show that, in order to maintain power as the
effect sizes shrink, the sample sizes need to grow at the rates claimed in Corollary 2.9.
To do this, we take the structural equation model parameterization from Richardson and
Spirtes (2002, Section 8); for Gk we set the weight of each edge on the discriminating path
to be ρs = 0.4 × 2−s , and of the other edges to be 0.5. We take a sample size of ninit × 22ks
for some ninit , so that this grows at the rate suggested by Corollary 2.9 to keep the power
constant in s. The plot in Figure 8 shows that the simulations agree with the predictions of
that result, since the power stays roughly constant as s grows.
Table 1 lists some of the colossal sample sizes that are needed to keep power constant in
longer discriminating paths as effect sizes shrink. In the k = 3 case, we have approximately
70% accuracy when ρ = 0.2 and n = 16 000, but to maintain this when ρ = 0.1 we need
n = 1 024 000. For k = 4 the respective sample sizes are an increase from n = 2 × 105 to
n = 5 × 107 (a 256-fold increase), and for k = 5 from n = 5 × 106 to n = 5 × 109 (a 1024-fold
increase).

5

Other Classes of Model

In this section we provide some further examples, involving both graphical and other models.

16

acc.
0.796
0.704
0.667
0.666
0.655
0.657
0.650

1.0
0.8
0.6

0.7

proportion correct

0.9

k=2
k=3
k=4
k=5

0

1

2

3

4

5

s

Figure 8: Ability to discriminate between Gk and Gk0 for various k and s. A sample of 2,500
data sets were drawn from the Gaussian graphical model according to the scheme described
in the text, with effect sizes ρs = 0.4 × 2−s ; the sample size was fixed at ninit × 22ks for some
ninit . The y-axis gives the proportion of times Gk correctly gave a lower deviance than Gk0 .
Each solid line corresponds to a value of k ∈ {2, 3, 4, 5} and a corresponding sample size
ninit ∈ {32, 250, 800, 5000} (dashed lines correspond to an initial sample size of 4ninit ). The
flattening out each line shows agreement with the prediction of Corollary 2.9. Note that the
vertical ordering of the different lines is merely a function of the choice of ninit . Larger values
of k and s were excluded because the relevant sample sizes grew too quickly (k = s = 5
corresponds to n = 5.6 × 1018 for the solid line).

17

5.1

Undirected Graphs

An undirected graphical model associates a simple undirected graph with a collection of
probability distributions. Under the pairwise Markov property, the model consists of distributions such that Xi ⊥
⊥ Xj | XV \{i,j} whenever i and j are not joined by an edge in the
graph. For Gaussian graphical models, this is equivalent to enforcing a zero in the (i, j)
entry in the inverse covariance matrix. It is a simple matter to see that no two distinct
undirected Gaussian graphical models ever overlap.
This helps to explain why undirected models are fundamentally easier to learn than
other classes, something which has been much exploited in high-dimensional statistics. For
example, the graphical lasso (Friedman et al., 2008; Witten et al., 2011) and neighbourhood
selection (Meinshausen and Bühlmann, 2006) methods allow very fast consistent model selection amongst undirected Gaussian graphical models, including in high-dimensional settings.
The result in Section 2.2 shows that a similar method could not exist for directed graphical
models without a prespecified ordering, at least not without failing in some regions of the
parameter space.

5.2

Time Series

The autoregressive moving average process of order (p, q), or ARMA(p, q) model, is a time
series model that assumes
Xt = εt +

p
X

φi Xt−i +

i=1

q
X

θi εt−i ,

t∈Z

i=1

where εt ∼ N (0, σ 2 ) and the parameters φ = (φi )pi=1 and θ = (θi )qi=1 are unknown. The
model is stationary and Gaussian, and therefore parameterized by σ 2 and the autocorrelations: i.e. γi ≡ Cor(Xt , Xt+i ) for i = 1, 2, . . .. Let the space spanned by γi be Li .
The special cases where p = 0 and q = 0 are respectively the MA(q) and AR(p) models.
These are identifiable, and the unique parameter values that lead to independence of the
Xt s are (φ, θ) = (0, 0). It is not hard to see that the derivative of the joint correlation
matrix with respect to either φi or θi is just Li . Hence the tangent space of an AR(p) or
MA(p) model at this point is of the form
i.i.d.

C(0,0) (Mp ) =

p
M

Li ,

i=1

and so the two AR(p) and MA(p) models are 1-equivalent at the point of joint independence. This suggests it will be hard to distinguish between these two types of process when
correlations are weak, though this may not matter if the aim of an analysis is predictive.

5.3

Nested Markov Models

Richardson et al. (2017) introduce nested Markov models, which are defined by a generalized
form of conditional independence that may hold under a Bayesian Network model with
hidden variables. For example, the two causal Bayesian networks in Figure 9 differ only in
the direction of a single edge (that between L and B). Assuming that the variable U is
unobserved, the model (a) implies a single observable conditional independence constraint:
Y ⊥
⊥ A | B.
The model in (b) does not imply any conditional independence constraint over the observed variables, but does imposes the restriction that
X
p(y | do(a, b)) ≡
p(l | a) · p(y | a, l, b)
l

18

U
A

L

B

U
Y

A

L

(a)

B

Y

(b)

Figure 9: Two causal models on five variables, differing in the direction of the L − B edge.
The variable U is assumed to be unobserved.

does not depend on a. This can be interpreted as the statement that A does not causally
affect Y , except through B. Note that, if L ⊥
⊥ B | A, then
X
p(y | do(a, b)) ≡
p(l | a) · p(y | a, l, b)
l

=

X

p(l | a, b) · p(y | a, l, b)

l

= p(y | a, b),
so the statement that this quantity does not depend upon a is the same as the ordinary
independence statement Y ⊥
⊥ A | B; a similar conclusion can be reached if Y ⊥
⊥ L | A, B.
Thus there are two distinct submodels along which the two models in Figure 9 intersect.
These two submodels also have linearly independent normal spaces because they correspond
to restrictions on p(a | y, b) and p(l | a, y, b). It follows that the two models are 2-nearequivalent at points where Y ⊥
⊥ A, L | B (where both submodels hold) by Theorem 3.2; note
that unlike for the case of ancestral graph models, this applies even in the case of discrete
variables.
So, even though in principle one can distinguish the two models in Figure 9 and determine
the orientation of the L − B edge, in practice the distinction may be hard to show with data.
This suggests it will be very difficult to learn the correct model without further information,
something borne out by the simulations in Shpitser et al. (2013). Interestingly, the approach
in that paper of setting higher-order parameters to zero in order to equalize parameter counts
between models may actually have made learning significantly harder, since in some cases
this means removing the only directions in the tangent space that differ between models.

6
6.1

Related Phenomena
Double Robustness

The phenomenon of double robustness of estimators has been exploited to improve estimation
in causal models (Scharfstein et al., 1999); see Kang and Schafer (2007) for an overview.
The essential element of double robustness
√ is to choose an estimating equation with two
components, such that the solution is a n-consistent estimator when one of the two parts
is correctly specified, even if the other is misspecified. For simplicity, we will just consider
examples where we assume independence rather than fitting a model.
Consider the simple causal model depicted in Figure 10, and suppose we are interested
in the causal effect of a binary variable X on the expectation of Y , but there is a continuous
P
measured confounder Z. The causal effect of X on Y is given by p(y | do(x)) P
= z p(z) ·
p(y | x, z), which is generally different from the ordinary conditional p(y | x) = z p(z | x) ·
p(y | x, z).
What happens if we use the ordinary conditional anyway? One can easily check that
p(y | x) = p(y | do(x)) if either X ⊥
⊥ Z or Y ⊥
⊥ Z | X, and hence the set of distributions
19

Z
X

Y

Figure 10: Simple causal model on three variables.

where our estimate is correct contains the union of points MX⊥
⊥Z ∪ MY ⊥
⊥Z|X . If we apply
Proposition 3.1 we find that any error in estimating the causal effect will be quadratic in
the distance from the point where Z ⊥
⊥ X, Y .
In this case the union does not arise as an intersection of two models, so we would not
view this as ‘equivalence’ of two models. However, if we instead consider the submodels given
by X ⊥
⊥ Y and X ⊥
⊥ Y [do(X)], their intersection is given by the intersection of MX⊥
⊥Y
with the union MX⊥
⊥Z ∪ MY ⊥
⊥Z|X ; these two models are therefore 2-near-equivalent.

6.2

Triple Robustness

Another phenomenon known as triple robustness 4 is observed in some causal models related
to mediation (Tchetgen and Shpitser, 2012). In this case, an estimator will be consistent
provided at least two out of three other quantities are correctly specified. Here we introduce
another result related to Proposition 3.1
Proposition 6.1. Suppose that f is a C 2 function and f (x) = 0 for any x ∈ Ni ⊂ M for
i = 1, . . . , m and that the tangent spaces Tθ (Ni ) jointly span M. Then, f (x) = O(kxk2 ).
Proof. In this case, for any differentiable function that is zero on M, whichever direction
we move in can be written as a linear combination of directions in the tangent spaces of
submodels Ni . It follows that the directional derivative of any such function is zero, in any
direction away from θ.
In light of this, suppose we have three submodels M1 , M2 , M3 each defined by constraints on linearly independent parts of the model M, and such that an estimator is consistent on the intersection of any two of them. The condition on the definition of the Mi means
that their normal spaces are linearly independent, so any vector is in the tangent space of
at least two such submodels. Their pairwise intersections thus satisfy the conditions of the
theorem, and the error in estimating the relevant parameter is quadratic in the distance to
the intersection M1 ∩ M2 ∩ M3 .

6.3

Post-double-Selection

Belloni et al. (2014) consider the problem of estimating a causal effect p(y | do(x)) in the
presence of a high-dimensional measured confounder ZI , I = {1, . . . , p} where p  n. We
can try to find a subset S ⊆ I such that S is much
P smaller than I, and ZS is sufficient to
control for the confounding, i.e. p(y | do(x)) = zS p(zS ) p(y | x, zS ). Formally, this will be
satisfied if we ignore any Zi such that Zi ⊥
⊥ X | Z−i or Zi ⊥
⊥ Y | X, Z−i .
In finite samples, selecting variables creates an omitted-variable bias, in which the decision
boundary of whether to drop a particular variable leads to a bias of order O(n−1/2 ) at some
points in the parameter space.
If we only exclude components of Zi for which both Zi ⊥
⊥ X | Z−i and Zi ⊥
⊥ Y | X, Z−i ,
then the order of the bias on our causal estimate is effectively squared and becomes O(n−1 );
this is because—for the same reason as in the discussion of double robustness—the bias
induced by a component Zi is at most quadratic in the distance of the true distribution
4 Perhaps

misleadingly, since it is strictly weaker than double robustness.

20

from the intersection of these two independences. Since this bias is small compared to the
sampling variance, it can effectively be ignored; this idea is referred by Belloni et al. (2014)
to as post-double-selection, and can be viewed as another consequence of the local geometry
of the model.

7

Algorithms for Learning Models with Overlap

Suppose we have a class of models Mi that overlap but are not 1-equivalent: that is, the
models all have different tangent cones. We have seen already that overlapping models place
restrictions on one class of computationally efficient methods, because they cannot be made
convex. This suggests that a method which attempts to learn the (linear) tangent cone
rather than selecting the model directly may be computationally advantageous. This can
be achieved by learning from a set of ‘surrogate’ models that have the same tangent spaces
as the original models, but that do not overlap.
To take a simple example, consider again the graphical models in Figure 2 for discrete
variables. The model M2 in (b) consists of the zero log-linear parameters: λXY = λXY Z =
0, while (a) involves log-linear parameters over the X, Y -margin M1 : λ0XY = 0. This
apparently makes model selection tricky because some models involve zeroes of ordinary
log-linear parameters, and some of marginal log-linear parameters.
However, one could try replacing M1 with a model that corresponds to the zero of the
ordinary log-linear parameter M01 : λXY = 0, as this has the same tangent space as M1
at the uniform distribution, and so it is ‘close’ to M1 in a precise sense. This suggests
that if we pursued a model selection strategy for ordinary log-linear parameters and learned
λXY = 0 but λXY Z 6= 0 (i.e. we chose M01 ), then we could conclude that this is sufficiently
similar to M1 to select this model from the graphical class.

7.1

Model Selection

Define Λi = hei i to be the vector space spanned by the ith coordinate axis. Suppose we
have a class of regular models Mi ⊆ Θ ⊆ Rk such that each model has a tangent space
at θ = 0 defined by a subset of coordinate axes: that is, for each model there is some set
c(Mi ) ⊆ {1, . . . , k} such that
M
C0 (Mi ) =
Λj .
j∈c(Mi )

Suppose further that our class of models is such that c(Mi ) 6= c(Mj ) for any i 6= j. We have
already seen that, if two models overlap, then it may be computationally difficult to learn
the correct one due to a lack of convexity. We will show that, under certain assumptions
about the true parameter being sufficiently close to θ = 0, we can learn the tangent space
itself, thereby circumventing the models’ lack of convexity.
Denote the sparsity pattern of a parameter by c(θ) = {i : θi 6= 0}; then θ ∈ T0 (M)
implies that c(θ) ⊆ c(T0 (M)) = c(M). However, note that c(θ) ∈ M does not imply this,
and in general the sparsity patterns of parameters in M is arbitrary.
Suppose we have a model selection procedure which returns a set Ŝ estimating S ⊆
{1, . . . , k} such that θS 6= 0 and θS c = 0. We will assume that the procedure is consistent, in
the sense that if θSnc = o(n−1/2 ) and |θsn | = ω(n−1/2 ) for each s ∈ S, we have P (Ŝ = S) → 1.
These conditions are satisfied by many common model selection methods such as BIC, or
an L1 -penalized selection method with appropriate penalty (and with Fisher information
matrix satisfying certain irrepresentability conditions, Nardi and Rinaldo, 2012).
The following result shows that we can adapt a model selection procedure of this kind
to learn models that overlap, by replacing each (possibly non-convex) model of interest by
a convex surrogate model with the same tangent space.

21

Theorem 7.1. Consider a sequence of parameters of the form θn = θ0 +n−γ h+O(n−2γ ) for
some 41 < γ < 12 , such that each θn ∈ M. Suppose our model selection procedure provides a
sequence of parameter estimates θ̃n .
Then
P (c(θ̃n ) = c(M)) → 1
as n → ∞.
Proof of Theorem 7.1. If i 6∈ c(M) then θni = O(n−2γ ) = o(n−1/2 ) since γ > 14 , while if
i ∈ c(M) then θni = Ω(n−γ ) = ω(n−1/2 ) since γ < 12 . By the conditions on our model
selection procedure then, we have the required consistency.
The condition γ > 14 is to ensure that the bias induced by using the surrogate model is
too small to detect at the specific sample size, and that our procedure will set the relevant
parameters to zero. Slower rates of convergence (i.e. 0 < γ ≤ 41 ) could also lead to a
satisfactory model selection procedure if we were simply to subsample our data or otherwise
‘pretend’ that n is smaller than it actually is. If γ > 12 , on the other hand, we will not have
asymptotic power to identify the truly non-zero parameters.
Note that P (c(θ̃n ) = c(θn )) does not tend to 1, since the sparsity pattern of the true
parameter is not the same as that of the tangent space of the model: it is merely ‘close’ to
having the correct sparsity.
The assumption that θn tends to 0 at the required rate may seem rather artificial:
some assumption of this form is unavoidable, simply because our results only hold in a
neighbourhood of points of intersection. The precise rate at which θn → 0 just needs to
be such that √
‘real’ effects do not disappear faster than we can statistically detect them (i.e.
slower than n), and that any other effects shrink fast enough that they are taken to be
zero. Any more ‘realistic’ framework would require conditions on the global geometry of the
models, and would be extremely challenging to verify.

7.2

Application to Bayesian Networks

Suppose that we have a sequence of discrete distributions pn , with kpn − p0 k = O(n−γ ) for
1
1
4 < γ < 2 , where p0 is the uniform distribution, and such that each pn is Markov with
respect to a Bayesian network G. Then a consistent method to determine G would be as
follows:
• using the log-linear lasso with penalty ν = nδ for
for pn ;

1
2

< δ < 1 to learn a log-linear model

• then, asymptotically, the sparsity pattern of the log-linear model is the same as that of
the original graph, so simply find the graph with skeleton given by i − j when λij 6= 0
and orient v-structures as i → k ← j if and only if λijk 6= 0.
This is far from an optimal approach, but does give an idea of how one might be able to
overcome the non-convexity inherent to overlapping models.

8

Discussion

We have proposed that the geometry of two models at points of intersection is a useful
measure of how statistically difficult it will be to distinguish between them, and shown that
when models’ tangent spaces are not closed under intersection this restricts the possibility of
using convex methods to perform model selection in the class. We have also given examples
of model classes in which this occurs and noted that in several cases, model selection is
indeed known to be difficult.
22

We suggest that special consideration should be given in model selection problems to
whether or not the class contains models that overlap or are 1-equivalent and—if it does—
to whether a smaller and simpler model class can be used instead; alternatively, additional
experiments may need to be performed to help distinguish between models. The results in
this paper provide a point of focus for new model selection methods and also for experimental
design. If we are able to work with a class of models that is less rich and therefore easier to
select from, then perhaps we ought to. If we cannot, it is useful to know in advance at what
points in the parameter space it is likely to be difficult to draw clear distinctions between
models, so that we can power our experiments correctly or just report that we do not know
which of several models is correct.
Acknowledgments
We thank Thomas Richardson for suggesting one of the examples, and several other readers
for helpful comments.

References
B. Aragam and Q. Zhou. Concave penalized estimation of sparse Gaussian Bayesian networks. Journal of Machine Learning Research, 16:2273–2328, 2015.
I. A. Beinlich, H. J. Suermondt, R. M. Chavez, and G. F. Cooper. The ALARM monitoring
system: A case study with two probabilistic inference techniques for belief networks. In
J. Hunter, J. Cookson, and J. Wyatt, editors, Proceedings of AIME 89: Second European Conference on Artificial Intelligence in Medicine., pages 247–256. Springer Berlin
Heidelberg, 1989.
A. Belloni, V. Chernozhukov, and C. Hansen. Inference on treatment effects after selection
among high-dimensional controls. The Review of Economic Studies, 81(2):608–650, 2014.
W. P. Bergsma and T. Rudas. Marginal models for categorical data. Annals of Statistics,
30(1):140–159, 2002.
C. M. Bishop. Pattern recognition and machine learning. Springer, 2007.
J. Bochnak, M. Coste, and M.-F. Roy. Real algebraic geometry, volume 36. Springer Science
& Business Media, 2013.
L. Breiman. Statistical modeling: The two cultures (with discussion). Statistical Science,
16(3):199–231, 2001.
P. Bühlmann and S. van de Geer. Statistics for High-Dimensional Data: Methods, Theory
and Applications. Springer, 2011.
P. Bühlmann, J. Peters, J. Ernest, et al. CAM: Causal additive models, high-dimensional
order search and penalized regression. Annals of Statistics, 42(6):2526–2556, 2014.
D. M. Chickering. Learning Bayesian networks is NP-complete. In Learning from data,
pages 121–130. Springer, 1996.
D. Cox, J. Little, and D. O’Shea. Ideals, Varieties, and Algorithms: An Introduction to
Computational Algebraic Geometry and Commutative Algebra. Springer, third edition,
2008.
D. R. Cox. Role of models in statistical analysis. Statistical Science, 5(2):169–174, 1990.
M. Drton. Discrete chain graph models. Bernoulli, 15(3):736–753, 2009.

23

M. Drton and S. Sullivant. Algebraic statistical models. Statistica Sinica, 17(4):1273–1297,
2007.
M. Drton, B. Sturmfels, and S. Sullivant. Lectures on algebraic statistics, volume 39. Springer
Science & Business Media, 2008.
M. Ferrarotti, E. Fortuna, and L. Wilson. Local approximation of semialgebraic sets. Annali
della Scuola Normale Superiore di Pisa, 1:1–11, 2002.
J. Friedman, T. Hastie, and R. Tibshirani. Sparse inverse covariance estimation with the
graphical lasso. Biostatistics, 9(3):432–441, 2008.
F. Fu and Q. Zhou. Learning sparse causal Gaussian networks with experimental intervention: regularization and coordinate descent. Journal of the American Statistical Association, 108(501):288–300, 2013.
J. Gu, F. Fu, and Q. Zhou. Adaptive penalized estimation of directed acyclic graphs from
categorical data. arXiv preprint arXiv:1403.2310, 2014.
J. D. Y. Kang and J. L. Schafer. Demystifying double robustness: A comparison of alternative strategies for estimating a population mean from incomplete data. Statistical Science,
22(4):523–539, 2007.
S. L. Lauritzen and N. Wermuth. Graphical models for associations between variables, some
of which are qualitative and some quantitative. Annals of Statistics, 17(1):31–57, 1989.
L. Le Cam. Locally asymptotically normal families of distributions. University of California
Publications in Statistics, 3:37–98, 1960.
N. Meinshausen and P. Bühlmann. High-dimensional graphs and variable selection with the
lasso. Annals of Statistics, 34(3):1436–1462, 2006.
Y. Nardi and A. Rinaldo. The log-linear group-lasso estimator and its asymptotic properties.
Bernoulli, 18(3):945–974, 2012.
Y. Ni, F. C. Stingo, and V. Baladandayuthapani. Bayesian nonlinear model selection for
gene regulatory networks. Biometrics, 71(3):585–595, 2015.
Pearl. Causality: models, reasoning and inference. CUP, second edition, 2009.
T. S. Richardson and P. Spirtes. Ancestral graph Markov models. Annals of Statistics, 30
(4):962–1030, 2002.
T. S. Richardson, R. J. Evans, J. M. Robins, and I. Shpitser. Nested Markov properties for
acyclic directed mixed graphs. arXiv:1701.06686, 2017.
G. V. Rocha, X. Wang, and B. Yu. Asymptotic distribution and sparsistency for `1 -penalized
parametric M-estimators with applications to linear SVM and logistic regression. arXiv
preprint arXiv:0908.1940, 2009.
T. Rudas, W. P. Bergsma, and R. Németh. Marginal log-linear parameterization of conditional independence models. Biometrika, 97:1006–1012, 2010.
D. O. Scharfstein, A. Rotnitzky, and J. M. Robins. Rejoinder to “adjusting for nonignorable
drop-out using semiparametric nonresponse models”. Journal of the American Statistical
Association, 94(448):1135–1146, 1999.
H. Schlichtkrull. Differentiable manifolds. 2011. Lecture Note of the Department of Mathematical Sciences, University of Copenhagen.

24

S. Shimizu, P. O. Hoyer, A. Hyvärinen, and A. Kerminen. A linear non-Gaussian acyclic
model for causal discovery. The Journal of Machine Learning Research, 7:2003–2030,
2006.
A. Shojaie and G. Michailidis. Penalized likelihood methods for estimation of sparse highdimensional directed acyclic graphs. Biometrika, 97(3):519–538, 2010.
I. Shpitser, R. Evans, T. Richardson, and J. Robins. Sparse nested Markov models with
log-linear parameters. In UAI-2013, pages 576–585, 2013.
A. Skrondal and S. Rabe-Hesketh. Generalized latent variable modeling: Multilevel, longitudinal, and structural equation models. Chapman & Hall/CRC Press, 2004.
P. Spirtes, C. Glymour, and R. Scheines. Causation, Prediction, and Search. MIT Press,
second edition, 2000.
E. J. T. Tchetgen and I. Shpitser. Semiparametric theory for causal mediation analysis:
efficiency bounds, multiple robustness, and sensitivity analysis. Annals of Statistics, 40
(3):1816, 2012.
R. Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society, Series B, 58:267–288, 1996.
C. Uhler, G. Raskutti, P. Bühlmann, and B. Yu. Geometry of the faithfulness assumption
in causal inference. Annals of Statistics, 41(2):436–463, 2013.
A. van der Vaart. The statistical work of lucien le cam. Annals of Statistics, 30(3):631–682,
2002.
A. W. van der Vaart. Asymptotic Statistics. Cambridge University Press, 1998.
D. M. Witten, J. H. Friedman, and N. Simon. New insights and faster computations for the
graphical lasso. Journal of Computational and Graphical Statistics, 20(4):892–900, 2011.
J. Zhang. On the completeness of orientation rules for causal discovery in the presence of
latent confounders and selection bias. Artificial Intelligence, 172(16):1873–1896, 2008.
P. Zwiernik, C. Uhler, and D. Richards. Maximum likelihood estimation for linear Gaussian
covariance models. Journal of the Royal Statistical Society, Series B, 2017. To appear.

A
A.1

Technical Results
Geometry

The proof below makes modest use of differential geometry; the basics may be found in
Schlichtkrull (2011).
Proof of Theorem 3.2. Let M1 be a C m manifold containing 0, and such that Ni are manifolds each with disjoint normal spaces. Let φ1 : U1 → Rd be a chart for M1 containing 0,
and such that the preimage of N1 ∩ M1 under φ1 consists of certain co-ordinate zeroes; e.g.
φ−1
1 (N1 ∩ M1 ) = {u ∈ U1 : u1 = · · · = uc1 = 0}.
Repeat for each Ni (using distinct coordinate zeros), and do the same for the corresponding
chart φ2 of M2 . Now, any component of either chart satisfies the properties of Proposition
3.1, and so we can write out its Taylor expansion to order m − 1 and the remainder is of
order O(kukm ).

25

Since whenever u1 = . . . = uc1 = 0 the models coincide, this shows that all derivatives of
the two charts are the same unless they involve at least one of these directions. This applies
to each of the other models in turn, so we find that derivatives of all orders are the same
apart from those containing at least one derivative from each model, which will be at least
an mth order partial derivative.
It follow that the distance between the two models close to 0 is of order O(εm ), and
hence they are m-near-equivalent.

A.2

Asymptotics

Definition A.1. Let (pθ : θ ∈ Θ) be a class of densities with respect to a measure µ indexed
by some open Θ ⊆ Rk . We say that this class is differentiable in quadratic mean (DQM) at
θ if there exists a vector `˙θ ∈ Rk such that
2
Z 
1 T˙ √
√
√
pθ+h − pθ − h `θ pθ dµ = o(khk2 ).
2
We further say that this class is uniformly differentiable in quadratic mean (UDQM) at θ if
there exists a neighbourhood N of θ such that
sup
θ∈N

Z 
√

pθ+h −

√

1
√
pθ − hT `˙θ pθ
2

2

dµ = o(khk2 ).

Note that
∂
`˙θ (x) =
log pθ (x),
∂θ
at such x that this is defined (which will be almost all x). Consequently it is almost surely
differentiable at a given θ, and the Fisher information I(θ) ≡ E`˙θ `˙Tθ is continuous in θ.
Formally uniform DQM seems to be strictly stronger than ordinary DQM, though we
are unable to construct any examples that satisfy the latter condition but not the former.
Such cases are likely to be rather pathological, and certainly most common classes (such as
exponential families) are UDQM at all θ.
Lemma A.2. Suppose (pθ : θ ∈ Θ) is UDQM at θ. Then for h̃ → 0,
Z h
i2
p
p
√
pθ+h − pθ+h̃ − (h − h̃)`˙θ+h̃ pθ+h̃ dµ = o(kh − h̃k2 ).
Proof. Note that the result is only meaningful if also h → 0. It is then just a direct
application of Definition A.1 with θ and h replaced by θ + h̃ and h − h̃ respectively.
Proposition A.3. Let Θ be an open√subset of Rk and pθ be a model which is uniformly
DQM at θ ∈ Θ. If kn = hn − h̃n and nkn → k with hn → 0, then
n

1
1 X T˙
k `θ+h̃n (Xi ) − k T I(θ)k + op (1)
`(θ + hn ) − `(θ + h̃n ) = √
2
n i=1


1 T
d
T
−→ N − k I(θ)k, k I(θ)k .
2
where I(θ) = E`˙θ `˙Tθ .
Proof. This is obtained by a simple adaptation of van der Vaart (1998, Theorem 7.2); we
make use of similar notation, so that Pθ̃n is abbreviated to P̃n , and Pθn and Pθ̃n have

26

˙ + h̃n ). First, we claim
densities pn and p̃n with respect to a measure µ. Let gn = knT `(θ
limn P̃n gn = 0; to see this, note
Z
p
1 p
P̃n gn =
gn p̃n 2 p̃n dµ,
2
√
√
√
1
and the fact that k pn − p̃n − 2 gn p̃n k2 = o(kkn k2 ) = o(n−1 ) (by Lemma A.2) means
that this has the same limit as
Z
p
p
√ √
√
n( pn − p̃n )( pn + p̃n ) dµ = 0.
Now define Wni = 2

hp

i
pn /p̃n (X ni ) − 1 , where X ni are independently generated from

P̃n . Then
Varθ̃n

n

n
X

1 X
gn (X ni )
Wni − √
n
i=1
i=1

!
≤ Eθ̃n
= 2n

√

2
nWni − gn (X ni )

Z 
√

pn −

p

p̃n − knT `˙θ+h̃n

p

p̃n

2

dµ

= 2no(kn2 ) = o(1),
by UDQM. Also,
Eθ̃n

n
X

!
Wni

Z
= 2n

√

pn


p
p̃n dµ − 1

i=1

= −n

Z 
√

pn −

p

p̃n

2

dµ

1
= − P̃n gn2 + o(1),
4
with the last equality again making use of Lemma A.2. Combining with the fact that
P̃n gn = 0, we find that the sequence of random variables
n
X

n

1 X
1
Wni − √
gn (X ni ) + P̃n gn2
4
n i=1
i=1

has mean and variance both tending to zero in probability; this implies that the random
variables themselves tend to zero in probability. The rest of the proof follows van der Vaart
(1998) exactly, and we obtain
n

1
1 X T˙
k `θ+h̃n (Xi ) + k T I(θ + h̃n )k + op (1).
`(θ + hn ) − `(θ + h̃n ) = √
2
n i=1
By continuity of the Fisher information, I(θ + h̃n ) = I(θ) + o(1), which gives the result.

A.3

Log-Linear Parameters

Let V be a finite set, and define the log-linear design matrix as a 2|V | × 2|V | matrix with
rows and columns indexed by subsets of V , such that
MA,B = (−1)|A∩B| .
We denote the Bth column (or equivalently row) of M by MB . Note that
K
MB =
M{v} ,
v∈B

27

where denotes the Hadamard (or point-wise) product. As an example, here is a log-linear
design matrix for three items.


1
1
1
1
1
1
1
1
∅
 1 −1
1 −1
1 −1
1 −1 
{1}


 1
1 −1 −1
1
1 −1 −1 
{2}


 {1, 2}
 1 −1 −1
1
1
−1
−1
1

M =
 1
{3}
1
1
1 −1 −1 −1 −1 



 1 −1
{1,
3}
1 −1 −1
1 −1
1 

{2,
3}

 1
1 −1 −1 −1 −1
1
1
{1,
2,
3}
1 −1 −1
1 −1
1
1 −1
For example, the fourth column of M is M{1,2} and is given by the pointwise product of the
second and third columns M{1} and M{2} . Note also that M is involutory—that is, its own
inverse—up to a constant: M −1 = 2−|V | M .
Let XV = (Xv )v∈V be a vector of binary random variables. We abbreviate the event
{Xv = 0 for all v ∈ C} to 0C , and similarly {Xv = 1 for all v ∈ C} to 1C . Let ηA =
log p(1A , 0V \A ). Then we define the log-linear parameters via the identities
λ = M −1 η,

η = M λ,

Letting p = (p(xV ) : xV ∈ XV ), assumed to be ordered in the same way as η so that
η = log p, we have
∂λ
∂λ ∂η
=
= M −1 (diag p(xV ))−1 .
∂p
∂η ∂p
Of interest to us is the connection between log-linear parameters within different marginal
distributions, known as marginal log-linear parameters (Bergsma and Rudas, 2002). Denote
the log-linear parameters within a marginal distribution XK by λK ≡ M K log p(xK ), where
M K is the appropriate restriction of M to rows and columns indexed by subsets of K. We
denote the ordinary log-linear parameter associated with a particular interaction set A by
λA = λVA .
Lemma A.4. The derivative of the parameter λK
A (with respect to p) lies in the span of the
⊥ XV \K . If K = A, then the converse also holds.
derivatives of λVA , . . . , λVV if XK ⊥
V
Additionally, if p(xV ) is uniform then λK
A and λA have the same derivative MA .
Proof. Let xV = (1B , 0V \B ). We have

and

∂λK
MA,B
A
=
∂p(xV )
p(xK )
X
∂λVAC
MA,B X
αC
=
αC MC,B
∂p(xV )
p(xV )

C⊆V \A

since MA∪C,B = MA,B · MC,B . For the derivatives
∂λV
AC
∂p

(3)

C⊆V \A

∂λK
A
∂p

to lie in the span of the derivatives

we need to find αC to solve
MA,B
p(xV )

X

αC MC,B =

C⊆V \A

MA,B
p(xK )

or equivalently
X

αC MC,B = p(xV \K | xK )

C⊆V \A

28

(4)

P
for each xV . If XV \K ⊥
⊥ XK then this becomes C⊆V \K αC MC,B = p(xV \K ), which has
a solution because it amounts to 2|V \K| linearly independent equations in 2|V \A| ≥ 2|V \K|
variables.
In the case that K = A, note that there are precisely as many variables as equations,
and since the coefficients MC,B expression given on the left of (4) do not vary with xK (since
C ∩ K = ∅), it is necessary for XV \K ⊥
⊥ XK in order for a solution to exist. For A ⊂ K a
similar argument shows that XV \K ⊥
⊥ XD | XK\D for some D ⊂ K with |D| ≤ |K \ A| is
sufficient.
If p(xV ) is uniform, then note that the derivative of λK
A does not depend upon K.
Since the map from λ to p is a smooth one, this yields us the following Corollary.
Corollary A.5. If ε ≡ kp − p0 k for some p0 under which all variables are independent, then
L
2
λK
A = λA + O(ε ).

A.4

Example: Discrete LWF Chain Graphs

Example A.6. Consider the graphs shown in Figure 11. Interpreted using the LauritzenFrydenberg-Wermuth (LWF) Markov property (Lauritzen and Wermuth, 1989), these three
graphs all represent distinct models. The graph in (a) satisfies the usual Markov property
for undirected graphs:
X1 ⊥
⊥ X4 | X2 , X3 ,

X2 ⊥
⊥ X3 | X1 , X4 .

The graphs in (b) and (c) satisfy these independences, as well as the respective independences
X1 ⊥
⊥ X2 (in the case of (b)), and X1 ⊥
⊥ X2 | X3 , X4 for (c). These two additional constraints
are generally distinct, but they coincide if any of the three remaining edges are not present:
i.e. if any of the conditional independences
X1 ⊥
⊥ X3 | X2 , X4 ,

X3 ⊥
⊥ X4 | X1 , X2 ,

X2 ⊥
⊥ X4 | X1 , X3 ,

also hold. Under either model, each of these constraints corresponds to a single zero loglinear parameter:
λ13 = 0,

λ34 = 0,

λ24 = 0.

Taking the models defined by these three constraints, we can apply Theorem 3.2 and find
that the two models in (b) and (c) are 3-near-equivalent at all points in the model of complete
independence.
One can extend this example arbitrarily by drawing a graph of the form 1 → 3 − 4 − · · · −
k ← 2 and comparing it to its undirected counterpart. In this case if the effect corresponding
to any of the k − 1 edges is missing, then the two models intersect. Hence, by Theorem 3.2
these two models are (k − 1)-near-equivalent.

29

1

2

1

2

1

2

3

4

3

4

3

4

(a)

(b)

(c)

Figure 11: Three chain graphs. When interpreted under the LWF Markov property, all
associated models satisfy X1 ⊥
⊥ X4 | X2 , X3 and X2 ⊥
⊥ X3 | X1 , X4 ; these constraints define
the model in (a). The submodel (b) additionally implies that X1 ⊥
⊥ X2 , whereas (c) implies
X1 ⊥
⊥ X2 | X3 , X4 .

30

