MAXIMUM LIKELIHOOD ESTIMATION IN GAUSSIAN MODELS
UNDER TOTAL POSITIVITY

arXiv:1702.04031v1 [stat.ME] 14 Feb 2017

STEFFEN LAURITZEN, CAROLINE UHLER, AND PIOTR ZWIERNIK

Abstract. We analyze the problem of maximum likelihood estimation for Gaussian
distributions that are multivariate totally positive of order two (MTP2 ). By exploiting
connections to phylogenetics and single-linkage clustering, we give a simple proof that
the maximum likelihood estimator (MLE) for such distributions exists based on n ≥ 2
observations, irrespective of the underlying dimension. Slawski and Hein [33], who
first proved this result, also provided empirical evidence showing that the MTP2
constraint serves as an implicit regularizer and leads to sparsity in the estimated
inverse covariance matrix, determining what we name the ML graph. We show that
the maximum weight spanning forest (MWSF) of the empirical correlation matrix is
a spanning forest of the ML graph. In addition, we show that we can find an upper
bound for the ML graph by adding edges to the MSWF corresponding to correlations
in excess of those explained by the forest. This also gives new theoretical results in
the study of inverse M-matrices. We provide globally convergent coordinate descent
algorithms for calculating the MLE under the MTP2 constraint which are structurally
similar to iterative proportional scaling. We conclude the paper with a discussion of
signed MTP2 distributions.

1. Introduction
Total positivity is a special form of positive dependence between random variables
that became an important concept in modern statistics; see, e.g., [3, 8, 20]. This
property (also called the MTP2 property) appeared in the study of stochastic orderings,
asymptotic statistics, and in statistical physics [14, 28]. Families of distributions with
this property lead to many computational advantages [2, 11, 30]. In a recent paper [12],
the MTP2 property was studied in the context of graphical models and conditional
independence in general. It was shown that MTP2 distributions have desirable Markov
properties. Our paper can be seen as a continuation of this work with a focus on
Gaussian distributions.
A p-variate real-valued distribution with density f w.r.t. a product measure µ is
multivariate totally positive of order 2 (MTP2 ) if the density satisfies
f (x)f (y) ≤ f (x ∧ y)f (x ∨ y).
A multivariate Gaussian distribution with mean µ and covariance matrix Σ is MTP2
if and only if the concentration matrix K := Σ−1 is an M-matrix, that is, Kij ≤ 0
Date: February 15, 2017.
2010 Mathematics Subject Classification. 60E15, 62H99, 15B48.
Key words and phrases. MTP2 distributions, Attractive Gaussian Markov random field (GMRF);
non-frustrated GRMF; Gaussian graphical model; inverse M-matrix; ultrametric.
1

2

MAXIMUM LIKELIHOOD IN TOTALLY POSITIVE GAUSSIAN MODELS

for all i 6= j or, equivalently, if all partial correlations are nonnegative. Such distributions were considered by Bølviken [5] and Karlin and Rinott [22]. Moreover, Gaussian
graphical models, or Gaussian Markov random fields, were studied in the context of
totally positive distributions in [26]. MTP2 Gaussian graphical models were shown to
form a sub-class of non-frustrated Gaussian graphical models, which themselves are a
sub-class of walk-summable Gaussian graphical models. Efficient structure estimation
algorithms for MTP2 Gaussian graphical models were given in [1] based on thresholding
covariances after conditioning on subsets of variables of limited size. Efficient learning
procedures based on convex optimization were suggested by Slawski and Hein [33] and
this paper is closely related to their approach; see also [4].
Throughout this paper, we assume that we are given n i.i.d. samples from N (µ, Σ),
where Σ is a positive definite matrix of size p × p. Without loss of generality, we assume
that µ = 0 and we focus on the estimation of Σ. We denote the sample covariance
matrix based on the n samples by S. Then the log-likelihood function is, up to additive
and multiplicative constants, given by
`(K; S) = log det K − tr(SK).
We denote the cone of real symmetric matrices of size p × p by Sp , its positive definite
elements by Sp0 , and its positive semidefinite elements by Sp0 . Note that `(K; S) is a
strictly concave function of K ∈ Sp0 . Since M-matrices form a convex subset of Sp0 ,
the optimization problem for computing the maximum likelihood estimator (MLE) for
MTP2 Gaussian models is a convex optimization problem. Slawski and Hein [33] showed
that the MLE exists, i.e., the global maximum of this optimization problem is attained,
when n ≥ 2. This yields a drastic reduction from n ≥ p without the MTP2 constraint.
In addition, they provided empirical evidence showing that the MTP2 constraint serves
as an implicit regularizer and leads to sparsity in the concentration matrix K.
In this paper, we analyze the sparsity pattern of the MLE K̂ under the MTP2
constraint. For a p × p matrix K we let G(K) denote the undirected graph on p nodes
with an edge uv if and only if Kuv 6= 0. Our main result is Theorem 4.2 that shows
that with probability one, the smallest maximum weight spanning forest (MWSF) of the
correlation matrix is a subgraph of the maximum likelihood (ML) graph G(K̂). This
provides a theoretical justification for initializing stepwise model selection methods
in the MWSF of the sample correlation matrix. In addition, in Proposition 4.4 we
obtain a simple upper bound for the ML graph G(K̂) by adding edges to the MWSF
corresponding to empirical correlations in excess of those provided by the MWSF. We
illustrate these results in the following example.
Example 1.1. We consider the carcass data that are discussed in [17] and can be
found in the R-library gRbase. This data set contains measurements of the thickness of
meat and fat layers at different locations on the back of a slaughter pig together with the
lean meat percentage on each of 344 carcasses. For our analysis we do not consider the
lean meat percentage, since, by definition, this variable should be negatively correlated
with the fat variables and positively correlated with the meat variables and hence the
joint distribution is unlikely to be MTP2 .

MAXIMUM LIKELIHOOD IN TOTALLY POSITIVE GAUSSIAN MODELS

3

The sample correlation matrix for this data set is given by
Fat11 Meat11 Fat12 Meat12 Fat13 Meat13

1.00
0.04
0.84
0.08 0.82 −0.03 Fat11
1.00
0.04
0.87 0.13
0.86Meat11
 0.04


0.84
0.04
1.00
0.01
0.83
−0.03

Fat12
S =

0.87
0.01
1.00 0.11
0.90Meat12
 0.08
 0.82
0.13
0.83
0.11 1.00
0.02Fat13
−0.03
0.86 −0.03
0.90 0.02
1.00 Meat13


and its inverse is
Fat11 Meat11 Fat12 Meat12 Fat13 Meat13

4.53
0.78 −2.34 −1.74 −1.72
1.01 Fat11
5.10 −0.25 −2.51 −0.76 −2.12Meat11
 0.78


−2.34 −0.25
4.52
1.40 −1.93 −0.94Fat12
=


1.40
7.00 −0.09 −4.11Meat12
−1.74 −2.51
−1.72 −0.76 −1.93 −0.09
4.11
0.53Fat13
1.01 −2.12 −0.94 −4.11
0.53
6.51 Meat13


S −1

This empirical covariance matrix is not MTP2 ; the positive entries in S −1 are highlighted in red. The MLE under MTP2 can be computed for example using cvx [16] in
matlab or using one of the simple coordinate descent algorithms discussed in Section 2.
The MLE of the correlation matrix, rounded to 2 decimals, is
Fat11 Meat11 Fat12 Meat12 Fat13 Meat13

1.00
0.10 0.84
0.09 0.82
0.09 Fat11
1.00 0.11
0.87 0.13
0.86Meat11
 0.10


0.84
0.11 1.00
0.09 0.83
0.09Fat12
Σ̂ =


0.87 0.09
1.00 0.11
0.90Meat12
 0.09
 0.82
0.13 0.83
0.11 1.00
0.11Fat13
0.09
0.86 0.09
0.90 0.11
1.00 Meat13


The entries of Σ̂ that changed compared to the sample correlation matrix S are highlighted in blue1. The sparsity pattern of K̂ = Σ̂−1 is captured by the ML graph G(K̂)
shown in Figure 1(a).
Note that all edges corresponding to blue entries in Σ̂ are missing in this graph. As
we show in Proposition 2.2, this is a consequence of the KKT conditions. Consider
now the maximum weight spanning forest of the complete graph with weights given
by the entries of S. In this example, the spanning forest is a chain represented by the
thick red edges in Figure 1(a). By Theorem 4.2 these edges form a spanning tree of the
ML graph G(K̂). by Proposition 4.4 the excess correlation graph (EC graph) yields
an upper bound for the ML graph G(K̂). In this example the EC graph identifies the
location of all non-edges apart from the edge Fat11 ∼ Meat12.
1We

note that Σ̂45 > S45 ; the entries appear equal only because of the 2-digit rounding.

4

MAXIMUM LIKELIHOOD IN TOTALLY POSITIVE GAUSSIAN MODELS

Meat11

Meat11

Fat13

Meat13

Fat12

Meat12

Fat11

(a) MTP2 constraint

Fat13

Meat13

Fat12

Meat12

Fat11

(b) Graphical lasso

Figure 1. Undirected Gaussian graphical models for carcass data. The
thick red edges in (a) correspond to the MWSF of the correlation matrix
and the blue edges in (b) to edges in the model chosen by the graphical
lasso which are not in the ML graph of the MTP2 solution.

Interestingly, applying various methods for model selection such as stepwise AIC,
BIC, or graphical lasso all yield similar graphs. In Figure 1(b) we show the graphical
model resulting from the graphical lasso, where the penalty parameter was chosen based
on a robustness analysis. Note that the graph obtained under the MTP2 constraint
as well as the upper bound are subgraphs of the model obtained using the graphical
lasso, indicating that ML estimation under the MTP2 constraint might be an interesting
alternative for model selection without the need of a tuning parameter. The extra edges
in Figure 1(b) are represented by thick blue lines.

The remainder of this paper is organized as follows: In Section 2, we review the duality theory that is known more generally for regular exponential families and specialize it
to MTP2 Gaussian distributions. This embeds the results by Slawski and Hein [33] into
the framework of exponential families and also leads to two related coordinate descent
algorithms for computing the MLE, one that acts on the entries of K and one that
acts on the entries of Σ. In Section 3, we show how the problem of ML estimation for
MTP2 Gaussian distributions is connected to single-linkage clustering and ultrametrics
as studied in phylogenetics. These observations result in a simple proof of the existence
of the MLE for n ≥ 2, a result that was first proven in [33]. Our proof is by constructing a primal and dual feasible point of the convex ML estimation problem for MTP2
Gaussian models. In Section 4, we prove our main results, namely the characterization
of a spanning forest of the ML graph G(K̂) and a simple upper bound for it. Finally,
in Section 5 we discuss how our results can be generalized to so-called signed MTP2
Gaussian distributions, where the distribution is MTP2 up to sign, i.e., |X| is MTP2 .
Such distributions were introduced by Karlin and Rinott in [21].

MAXIMUM LIKELIHOOD IN TOTALLY POSITIVE GAUSSIAN MODELS

5

2. Duality theory for maximum likelihood estimation under MTP2
We start this section by formally introducing MTP2 distributions and then discuss
the duality theory for Gaussian MTP2 distributions. Let V := {1, 2, . . . , p} be a finite
set and
Q let X = (Xv , v ∈ V ) be a random vector with density f on the product space
X = v∈V Xv , where Xv ⊆ R is the state space of Xv , inheriting the order from R. We
define the coordinate-wise minimum and maximum as
x ∧ y = (min(xv , yv ), v ∈ V ),

x ∨ y = (max(xv , yv ), v ∈ V ).

Then we say that X or the distribution of X is multivariate totally positive of order
two (MTP2 ) if its density function f on X satisfies
(1)

f (x)f (y)

≤

f (x ∧ y)f (x ∨ y)

for all x, y ∈ X .

In this paper, we concentrate on the Gaussian setting. It is easy to show that a
Gaussian distribution with mean µ and covariance matrix Σ is MTP2 if and only if
K = Σ−1 is an M-matrix, i.e.,
(i) Kvv > 0 for all v ∈ V ,
(ii) Kuv ≤ 0 for all u, v ∈ V with u 6= v.
Properties of M-matrices were studied by Ostrowski [29] who chose the name to honor
H. Minkowski. The connection to multivariate Gaussian distributions was established
by Bølviken [5] and by Karlin and Rinott [22].
We denote the set of all M-matrices of size p × p by Mp . Note that Mp is a convex
cone. In fact, it is obtained by intersecting the positive definite cone Sp0 with all the
coordinate half-spaces
p
Hij
= {X ∈ Sp | Xij ≤ 0}
with i 6= j. For a convex cone C we denote its closure by C. Then Mp is given
p
and the ML estimation problem for Gaussian MTP2 models can be
by Sp0 ∩i<j Hij
formulated as the following optimization problem:
maximize
(2)

K

log det(K) − trace(KS)

subject to K ∈ Mp

This is a convex optimization problem, since the objective function is concave on Sp0 .
Next, we introduce a second convex cone N p that plays an important role for ML
estimation in Gaussian MTP2 models. To formally define this cone, we introduce two
partial orders on matrices. Let A, B be two p × p matrices. Then A ≥ B means that
Aij ≥ Bij for all (i, j) ∈ V × V , and A  B means that A − B ∈ Sp0 . Then the cone
N p is defined as the negative closure of Sp0 , i.e.
N p = {X ∈ Sp | ∃Y ∈ Sp0 with X ≤ Y and diag(X) = diag(Y )}.
To simplify notation, we will suppress the dependence on p and write S, S0 , S0 , M
and N , when the dimension is clear. In the following result, we show that the cones N
and M are dual to each other.

6

MAXIMUM LIKELIHOOD IN TOTALLY POSITIVE GAUSSIAN MODELS

Lemma 2.1. The closure of N is the dual to the cone of M-matrices M, i.e.

(3)
N = S ∈ S | hS, Ki ≥ 0 for all K ∈ M .
Proof. We denote the dual of a convex cone C by C ∨ . Let C1 , C2 be two convex cones.
Then it is an easy exercise to verify that
(C1 ∩ C2 )∨ = C1∨ + C2∨ ;

(4)

here + denotes the Minkowski sum. Note that
S∨0 = S0

and

∨
Hij
= Hij .

This completes the proof, since M = S0 ∩i<j Hij and (4) can be applied inductively
to any finite collection of convex cones.

Using the cones M and N we now determine conditions for existence of the MLE in
Gaussian MTP2 models and give a characterization of the MLE. We say that the MLE
does not exist if the likelihood does not attain the global maximum.
Proposition 2.2. Consider a Gaussian MTP2 model. Then the MLE Σ̂ (and K̂)
exists for a given sample covariance matrix S on V if and only if S ∈ N . It is then
equal to the unique element Σ̂  0 that satisfies the following system of equations and
inequalities
(5)

(Σ̂−1 )uv ≤ 0 for all u 6= v,

(6)

Σ̂vv − Svv = 0 for all v ∈ V ,

(7)

(Σ̂uv − Suv ) ≥ 0 for all u 6= v,

(8)

(Σ̂uv − Suv )(Σ̂−1 )uv = 0 for all u 6= v,

Proof. The Lagrangian to (2) is given by
L(K, λ) = log det(K) − trace(KS) −

X

λuv Kuv ,

u<v

where we omitted the constraint K ∈ S0 . This can be done, since it is assumed that
K is in the domain of L. Maximizing L with respect to K we obtain

Svv
for u = v ∈ V
−1
(K̂ )uv =
Suv + λuv otherwise.
The Lagrange dual function is obtained by plugging in K̂ for K in L, i.e.,
g(λ) = log det(K̂) − tr(K̂ −1 K̂) = log det(K̂) − p.
Denoting K̂ −1 by Σ̂ and since λij ≥ 0, the dual optimization problem to (2) is given by
minimize
Σ0

(9)

− log det(Σ) − p

subject to Σvv = Svv , for all v ∈ V ,
Σuv ≥ Suv , for all u 6= v,

Note that the identity matrix is a strictly feasible point for (2). Hence, the MLE
does not exist if and only if the likelihood is unbounded. Since by Slater’s constraint

MAXIMUM LIKELIHOOD IN TOTALLY POSITIVE GAUSSIAN MODELS

7

qualification strong duality holds for the optimization problems (2) and (9), the MLE
does not exist if and only if S ∈
/ N.
Next we characterize the MLE. By the KKT conditions for convex optimization
problems, the primal-dual optimal points are characterized by
(1)
(2)
(3)
(4)

primal feasibility resulting in (5),
dual feasibility resulting in (6) and (7),
complimentary slackness resulting in (8),
stationarity of the Lagrangian resulting in tr(S Σ̂−1 ) = p, which is already implied by (8).

This completes the proof.



Remark 2.3. Proposition 2.2 can easily be extended to provide properties for the existence of the MLE and a characterization of the MLE for Gaussian graphical models
under MTP2 . In this case, let G = (V, E) be an undirected graph. Then the primal
problem has additional equality constraints, namely Kuv = 0 for all uv ∈
/ E, and hence
the inequality constraints in the dual problem are restricted to the entries in E, i.e.,
Σuv ≥ Suv for all uv ∈ E. Note that if the MLE of Σ based on S exists in the Gaussian
graphical model over G, it also exists in the Gaussian graphical model over G under
MTP2 , since without the MTP2 constraint the MLE needs to satisfy Σ̂uv = Suv for
all uv ∈ E. In fact, as shown by Slawski and Hein [33], see also Theorem 3.5 below,
the MLE under MTP2 exists for all sample sizes n ≥ 2 and this clearly remains so if
graphical model constraints are added.

We define the maximum likelihood graph (ML graph) Ĝ to be the graph determined
by K̂, i.e. Ĝ = G(K̂), where K̂ = Σ̂−1 is the MLE of K under MTP2 . We then have
the following important corollary of Proposition 2.2.
Corollary 2.4. Consider the Gaussian graphical model determined by Kuv = 0 for
uv 6∈ E(Ĝ), where Ĝ is the ML graph under MTP2 . Let Σ̄ be the MLE of Σ under that
Gaussian graphical model (without the MTP2 constraint). Then Σ̂ = Σ̄.
Proof. The MLE of Σ under the Gaussian graphical model with graph Ĝ is the unique
element Σ̄  0 that satisfies the following system of equations:
Σ̄vv − Svv = 0 for all v ∈ V ,
Σ̄uv − Suv = 0 for all uv ∈ E(Ĝ),
(Σ̄−1 )uv = 0 for all uv 6∈ E(Ĝ).
Proposition 2.2 says that also Σ̂ satisfies these equations and hence we must have
Σ̄ = Σ̂.

Note that this corollary highlights the role of the complementary slackness condition
(8) in inducing sparsity of the MTP2 solution.
We emphasize that the MLE under MTP2 is equivariant w.r.t. changes of scale so
that without loss of generality we can assume that the sample covariance is normalized,

8

MAXIMUM LIKELIHOOD IN TOTALLY POSITIVE GAUSSIAN MODELS

i.e. Sii = 1 or, equivalently, S = R, where R is the correlation matrix. For certain of
the subsequent developments this represents a convenient simplification.
Lemma 2.5. Let S be the sample covariance matrix, R the corresponding sample correlation matrix. Denote by Σ̂S and Σ̂R the MLE in Proposition 2.2 based on S and R,
respectively. Then
p
Σ̂Sij =
for all i, j ∈ V.
Sii Sjj Σ̂R
ij
Proof. Denote by D a diagonal matrix such that Dii =
The likelihood function based on S is

√

Sii , so that S = DRD.

log det K − tr(SK) = log det K − tr(RDKD).
P
If K 0 = DKD, this can be rewritten as log det K 0 − tr(RK 0 ) − i log Sii . Therefore, if K̂ R is the maximizer of log det K − tr(RK) under the MTP2 constraints, then
D−1 K̂ R D−1 is also an M-matrix and the maximizer of log det K − tr(SK).

We end this section by providing simple coordinate descent algorithms for ML estimation under MTP2 . Although interior point methods run in polynomial time, for
very large Gaussian graphical models it is usually more practical to apply coordinate
descent algorithms. In Algorithms 1 and 2 we describe two methods for computing the
MLE that only use optimization problems of size 2 × 2 which have a simple and explicit
solution, and iteratively update the entries of K, respectively of Σ. Algorithms 1 and 2
are inspired by the corresponding algorithms for Gaussian graphical models; see, for
example, [10, 35, 37]. Slawski and Hein [33] also provide a coordinate descent algorithm
for estimating covariance matrices under MTP2 . However, their method updates one
column/row of Σ at a time.
We first analyze Algorithm 1. Let A = {u, v} and B = V \ A. Then note that
the objective function can be written in terms of the 2 × 2 Schur complement K 0 =
−1
KAA − KAB KBB
KBA , since up to an additive constant
log det K − trace(KS) = log det K 0 − trace(K 0 SAA ).
Algorithm 1 Coordinate descent on K.
Input: Sample covariance matrix S, and precision .
Output: MLE K̂ ∈ M.
1. Let K 0 := K 1 := (diag(S))−1 .
2. Cycle through entries u 6= v and solve the following optimization problem:
maximize log det(K) − trace(KS)
K0

subject to Kuv ≤ 0,
Kij = Kij1 for all ij ∈ (V × V ) \ {uu, vv, uv},
and update K 1 = K.
3. If ||K 0 − K 1 ||1 < , set K̂ = K 1 . Otherwise, set K 0 = K 1 and return to 2.

MAXIMUM LIKELIHOOD IN TOTALLY POSITIVE GAUSSIAN MODELS

9

Algorithm 2 Coordinate descent on Σ.
Input: Sample covariance matrix S  0, and precision .
Output: MLE Σ̂ with Σ̂−1 ∈ M.
1. Let Σ0 := Σ1 := S
2. Cycle through entries u 6= v and solve the following optimization problem:
maximize log det(Σ)
Σ0

subject to Σuv ≥ Suv ,
Σij = Σ1ij for all ij ∈ (V × V ) \ {uv}.
and update Σ1 = Σ.
3. If ||Σ0 − Σ1 ||1 < , set Σ̂ = Σ1 . Otherwise, set Σ0 = Σ1 and return to 2.
−1
Defining L := KAB KBB
KBA , then the optimization problem in step (2) of Algorithm 1
is equivalent to
maximize
log det(K 0 ) − trace(K 0 SAA )
0
K 0

0
subject to K12
+ L12 ≤ 0.
−1
The unconstrained optimum to this problem is given by K 0 = SAA
and is attained if
−1
and only if (SAA )12 + L12 ≤ 0, or equivalently, if and only if

L12 ≤

Suv
.
2
Suu Svv − Suv

0
Otherwise the KKT conditions give that K12
= −L12 . Maximizing over the remaining
0
two entries of K leads to a quadratic equation, which has one feasible solution
p
p
2
1
+
4S
S
L
1
+
1 + 4Suu Svv L212
1
+
uu
vv
12
0
0
0
(10) K11
, K22
=
, K12
= −L12 .
=
2Suu
2Svv

Then the solution to the optimization problem in step (2) is given by KAA = K 0 + L.
Dual to this algorithm, one can define an algorithm that iteratively updates the offdiagonal entries of Σ by maximizing the log-likelihood in direction Σuv and keeping
all other entries fixed. This procedure is shown in Algorithm 2. Similarly as for
Algorithm 1, the solution to the optimization problem in step (2) can be given in
closed-form. Defining A = {u, v}, B = V \ A and L = ΣAB Σ−1
BB ΣBA , then analogously
as in the derivation above, one can show that the solution to the optimization problem
in step (2) of Algorithm 2 is given by
(11)

Σuv = max{Suv , L12 }.

We end by proving that Algorithms 1 and 2 indeed converge to the MLE. We here
assume that n ≥ 2 to guarantee existence of the MLE. Note that the suggested starting
points for both algorithms can be modified.
Proposition 2.6. Algorithms 1 and 2 converge to the MLE K̂ = Σ̂−1 ∈ M.

10

MAXIMUM LIKELIHOOD IN TOTALLY POSITIVE GAUSSIAN MODELS

Proof. The convergence to the MLE is immediate for Algorithm 2 because it is a
coordinate descent method applied to a smooth and strictly concave function; see,
e.g., [25]. For Algorithm 1 we use the fact that it is an example of iterative partial
maximization. To prove convergence to the MLE we we will show that the assumptions
of Proposition A.3 in [23] hold. The log-likelihood function that we are trying to
maximize is strictly concave and so the maximum is unique. Clearly K is the maximum
if and only if it is a fixed point of each update. It only remains to show that updates
depend continuously on the previous value. For a given S fix K and consider a sequence
of points Kn converging to K. Denote by K̃ and K̃n the corresponding one-step updates.
We want to show that K̃n also converges to K̃. As above, let A = {u, v}, B = V \ A,
−1
−1
K 0 = KAA − KAB KBB
KBA and L = KAB KBB
KBA . Outside of the block K̃AA this
convergence is trivial; so we focus only on the three entries in K̃AA . The function
0
0
0
L12 7→ (K11
, K22
, K12
) is continuous if and only if each coordinate is. It is clear that
these functions are continuous if L12 6= Suu SSvvuv−S 2 . It remains to show that if L12 =
uv
−1
Suv
0
the
update
in
(10)
gives
K
=
S
,
which
can be easily checked directly. 
2
AA
Suu Svv −S
uv

3. Ultrametric matrices and inverse M-matrices
In this section we exploit the link to ultrametrics in order to construct an explicit
primal and dual feasible point of the maximum likelihood estimation problem. A nonnegative symmetric matrix U is said to be ultrametric if
(i) Uii ≥ Uij for all i, j ∈ {1, . . . , p},
(ii) Uij ≥ min{Uik , Ujk } for all 1 ≤ i, j, k ≤ p.
The connection to M-matrices is established by the following result; see [9, Theorem 3.5].
Theorem 3.1. Let U be an ultrametric matrix with strictly positive entries on the
diagonal. Then U is nonsingular if and only if no two rows are equal. Moreover, if U
is nonsingular then U −1 is an M-matrix.
The main reason why ultrametric matrices are relevant in this paper is the following
construction, which is similar to constructions used in in phylogenetics [31, Section 7.2]
and single linkage clustering [15].
Let R be a symmetric p × p positive semidefinite matrix such that Rii = 1 for all
i = 1, . . . , p. Consider the graph G+ over V with an edge between i and j whenever
Rij is positive. Assign to each edge the corresponding positive weight Rij and note
that G+ in general does not have to be connected. Define a p × p matrix Z by setting
Zii = 1 for all i and
(12)

Zij := max min Ruv ,
P

uv∈P

where the maximum is taken over all paths in G+ between i and j and is set to zero if
no such path exists. We call Z the single-linkage matrix based on R.

MAXIMUM LIKELIHOOD IN TOTALLY POSITIVE GAUSSIAN MODELS

11

Example 3.2. Suppose that



1
−0.5 0.5 0.6
−0.5
1
0.4 −0.1

R = 
 0.5
0.4
1
0.2 
0.6 −0.1 0.2
1
Then G+ and Z are given by

1
0.6

2

4

0.5
0.4


1 0.4 0.5 0.6
0.4 1 0.4 0.4

Z = 
0.5 0.4 1 0.5 .
0.6 0.4 0.5 1


0.2

3
For example, to get Z12 we consider two paths 1 − 3 − 2 and 1 − 4 − 3 − 2. The
minimum of Ruv over the first path is 0.4 and over the second path 0.2. This gives
Z12 = 0.4. Note that Z ≥ R, Z is invertible, and Z −1 is an M-matrix. We now show
that this is an example of a more general phenomenon.

Proposition 3.3. Let R be a symmetric p × p positive semidefinite matrix satisfying
Rii = 1 for all i = 1, . . . , p. Then the single-linkage matrix Z based on R is an
ultrametric matrix with Zij ≥ Rij for all i 6= j. If, in addition, Rij < 1 for all i 6= j,
then Z is nonsingular and therefore an inverse M-matrix.
Proof. We first show that Z is an ultrametric matrix. Z is symmetric by definition.
Because R is positive semidefinite, Rij ≤ 1 for all i, j and from (12) it immediately
follows that Zij ≤ 1 and therefore Zii ≥ Zij for all i, j as needed. Finally, to prove
condition (ii) in the definition of ultrametric, let i, j, k ∈ V . Suppose first that i, j, k
lie in the same connected component of G+ . Let P1 , P2 be the paths in G+ such that
Zik = minuv∈P1 Ruv and Zjk = minuv∈P2 Ruv . Let P12 be the path between i and j
obtained by concatenating P1 and P2 . Then
Zij = max min Ruv ≥
P

uv∈P

min Ruv = min{Zik , Zjk }.

uv∈P12

Now suppose that i, j, k are not in the same connected component of G+ . In that case
0 ∈ {Zij , Zik , Zjk }. Because zero is attained at least twice, again Zij ≥ min{Zik , Zjk }.
Hence, Z is an ultrametric matrix. The fact that Zij ≥ Rij for all i, j follows directly
by noting that the edge ij forms a path between i and j.
Suppose now that Rij < 1 for all i 6= j. In that case also Zij < 1 for all i 6= j.
From this it immediately follows that no two rows of Z can be equal. Indeed, if the
i-th row is equal to the j-th row for some i 6= j, then necessarily Zij = Zii = Zjj , a
contradiction. From Theorem 3.1 it follows that Z is an inverse M-matrix.

As a direct consequence we obtain the following result.

12

MAXIMUM LIKELIHOOD IN TOTALLY POSITIVE GAUSSIAN MODELS

Proposition 3.4. Let S be a symmetric positivep
semidefinite matrix with strictly positive entries on the diagonal and such that Sij < Sii Sjj . Then there exists an inverse
M-matrix Z such that Z ≥ S and Zii = Sii for all i = 1, . . . , p.
Proof.
p Apply Proposition 3.3 to the normalized version R of S, with entries Rij :=
Sij / Sii Sjj . Because Rij < 1 for all i 6= j, the corresponding single-linkage matrix Z 0 is
p
ultrametric with Z 0 ≥ R and Z 0 is an inverse M-matrix. Define Z by Zij = Sii Sjj Z 0 .
Then Z ≥ S and Zii = Sii for all i = 1, . . . , p. Moreover, Z is an inverse M-matrix
because Z 0 is.

Proposition 3.4 is very important for our considerations. A basic application is an
elegant alternative proof of the main result of [33], which says that the MLE under
MTP2 exists with probability one as long as n ≥ 2. This is in high contrast with
the existence of the MLE in Gaussian graphical models without additional constraints;
see [36].
Theorem 3.5 (Slawski and Hein [33]). p
Consider a Gaussian MTP2 model and let S be
the sample covariance matrix. If Sij < Sii Sjj for all i 6= j then the MLE Σ̂ (and K̂)
exists and it is unique. In particular, if the number n of observations satisfies n ≥ 2,
then the MLE exists with probability 1.
Proof. The sample covariance matrix is a positive semidefinite matrix with strictly
positive diagonal entries. We can apply Proposition 3.4 to obtain an inverse M-matrix
Z that satisfies Z ≥ S and Zii = Sii for all i. It follows that Z satisfies primal feasibility
(5) and dual feasibility (6) and (7). By Proposition 2.2 the MLE exists and it is unique
by convexity of the problem.

The matrix Z can be computed efficiently2. To see that, note first that in Example 3.2
we could first consider the chain T of the form 2 − 3 − 1 − 4, which is the maximal
weight spanning forest of G+ and then construct Z by
Zij = min Ruv ,
uv=ij

where ij denotes the unique path between i and j in T . For example Z12 = 0.4, which
corresponds to the minimal weight on the path 2−3−1. This is a general phenomenon.
Suppose again that R is a symmetric p × p positive semidefinite matrix satisfying
Rii = 1 for all i = 1, . . . , p. Let MWSF(R) be the set of all minimal maximum weight
spanning forests of R. Note that all edge weights of any such forest F ∈ MWSF(R)
must be positive so we must have F ⊆ G+ . Also, if R is an empirical correlation
matrix, then MWSF(R) will be a singleton with probability one and in such cases we
shall mostly speak of the MWSF.
Proposition 3.6. The single-linkage matrix Z as defined in (12) is block diagonal with
blocks corresponding to the connected components of any F ∈ MWSF(R). Within each
block all elements are strictly positive and given by
Zij = min Ruv ,
uv∈ij
2In

our computations we use the single-linkage clustering method in R.

MAXIMUM LIKELIHOOD IN TOTALLY POSITIVE GAUSSIAN MODELS

13

where ij is the unique path between i and j in a maximal weight spanning tree of R. In
particular, Zij = Rij for all edges of MWSF(R).
Proof. First suppose that i, j ∈ V lie in two different components of F ∈ MWSF(R).
This means that there is no path between i and j in G+ and so, by definition, Zij = 0.
Because Zij > 0 if i, j lie in the same component of F , Z is block diagonal with blocks
corresponding to connected components of MWSF(R).
The rest of the proof is an adaptation of a proof of a related result [31, Proposition
7.2.10]. Suppose that i, j ∈ V lie in the same connected component of F and denote
the tree in F corresponding to this component by T . By definition Zij ≥ minuv∈ij Ruv .
Suppose that Zij > minuv∈ij Ruv . We obtain the contradiction by showing that under
this assumption T cannot be a maximum weight spanning tree of the corresponding
connected component of G+ . Let kl be a minimum weight edge in the unique path
between i and j in T . Since Zij > Rkl , there exists a path P in G+ between i and
j such that Ruv > Rkl for every uv in P . Now deleting kl from T partitions the
corresponding connected component of G+ into two sets with i being in one and j
being in the other block. Since P connects i and j in G+ , there must be an edge k 0 l0
(distinct from kl) in P whose end vertices lie in different blocks of this partition. Let T 0
be the spanning tree obtained from T by deleting kl and adding k 0 l0 . Since Rk0 l0 > Rkl ,
the total weight of T 0 is greater than T , which is a contradiction. We conclude that
Zij = minuv∈ij Ruv for all i, j in the same connected component of G+ .

To conclude this section, we note that the starting point Σ0 of Algorithm 2 is arbitrary
as long as Σ0  0 and Σ0 ≥ S. The single-linkage matrix constitutes another generic
choice when S = R is used as input. This is a particularly desirable starting point,
since it can also be used when p > n, in which case R ∈
/ S0 and hence not feasible as
a starting point.
4. The maximum likelihood graph
Fitting a Gaussian model with MTP2 constraints tends to induce sparsity in the
maximum likelihood estimate K̂. In this section, we analyze the sparsity pattern that
arises in this way. We assume again without loss of generality that S = R is a sample
correlation matrix so that Rii = 1 for all i and Rij < 1 for all i 6= j. Consider again
the weighted graph G+ on p nodes with edges corresponding precisely to the positive
off-diagonal entries of R with weights Rij . We begin this section with a basic lemma
that reduces our analysis to the case where the graph G+ is connected.
Lemma 4.1. The MLE Σ̂ under MTP2 is a block diagonal matrix with strictly positive
entries in each block. The blocks correspond precisely to trees in MWSF(R).
Proof. Firstly, since Σ̂ is an inverse M-matrix, it is block diagonal with strictly positive entries in each block; see, e.g., Theorem 4.8 in [19]. We will show that each block
of Σ̂ corresponds precisely to a tree in MWSF(R).
Denote the vertex sets for a forest F ∈ MWSF(R) as T1 , . . . , Tk and the blocks of Σ̂
as B1 , . . . , Bl . Firstly, for any Ti there must be a j so that Ti ⊆ Bj ; this is true since

14

MAXIMUM LIKELIHOOD IN TOTALLY POSITIVE GAUSSIAN MODELS

all entries in R along the edges of Ti are positive and thus Σ̂uv ≥ Ruv > 0. Thus the
block partitioning corresponding to the trees is necessarily finer than that of Σ̂.
On the other hand, suppose that two different trees Ti and Tj in F are in the same
block of Σ̂ so that Σ̂uv > 0 for all u ∈ Ti and v ∈ Tj . Then, as we must have Ruv ≤ 0,
also necessarily Σ̂uv − Ruv > 0. Complementary slackness (8) now implies that K̂uv = 0
for all u ∈ Ti and v ∈ Tj , and hence K̂ is block-diagonal with blocks corresponding to
the trees in F . Since Σ̂ = K̂ −1 , we also get Σ̂uv = 0 which contradicts that u and v are
in the same block of Σ̂.

This result shows that, without loss of generality, we can always assume that G+
is connected and then MWSF(R) = MWST(R) consists of trees only. If there are
more than one connected component, we simply compute the MLE for each component
separately and combine them together in block diagonal form. Hence, from now on we
always assume that all forests in MWSF(R) are just trees.
The following two theorems are the main results of the present paper.
Theorem 4.2. Let Σ̂ be the MLE of the correlation matrix under MTP2 based on the
sample correlation matrix R. Then any F ∈ MWSF(R) is a subgraph of the ML graph
Ĝ = G(K̂) with probability one.
Proof. Let Σ̄ denote a solution to the following optimization problem:
maximize log det(Σ)
Σ

subject to Σii = 1, for all i = 1, . . . , p,
Σij ≥ Rij , for all i 6= j,
Σij = Rij , for all ij ∈ MWSF(R),
Σ is an inverse M-matrix,
Σ  0.
Note that Σ̄ exists since the single linkage matrix Z is a feasible solution according to
Proposition 3.3 and Proposition 3.6, and the objective function is bounded because the
maximum is attained without the constraint that Σ is an inverse M-matrix.
Next we want to show that Σ̂ = Σ̄. We prove this by applying a step of Algorithm 2 on
an entry kl ∈ F ∈ MWSF(R). Without loss of generality we assume that kl = 12 ∈ F
and fix B = V \{1, 2}. Because Σ̄ is an inverse M-matrix, Σ̄12 ≥ Σ̄1B (Σ̄BB )−1 Σ̄B2 = L12 .
Also Σ̄12 = R12 by construction and so R12 ≥ L12 . It follows that the update (11) makes
no change in Σ̄, which proves that Σ̂ = Σ̄.
The last argument shows that Σ̂ij = Rij for any edge in a maximal weight spanning
forest of R. In other words, by Corollary 2.4, Σ̂ = Σ̄ is also the MLE under the
Gaussian graphical model with graph Ḡ = MWSF(R) ∪ G(K̂). However, this happens
with probability zero unless MWSF(R) ⊆ G(K̂) and hence the conclusion follows. 
At the heart of the proof of Theorem 4.2 we showed the following remarkable property
of the dual convex optimization problems (2) and (9): dual constraints (i.e., Σij ≥ Rij )

MAXIMUM LIKELIHOOD IN TOTALLY POSITIVE GAUSSIAN MODELS

15

that are active at an initial point that is primal feasible (i.e., Σ is an inverse M-matrix)
stay active at the solution. The single linkage matrix Z is an example of an inverse Mmatrix satisfying Zij = Rij for all ij ∈ MWSF(R) and hence MWSF(R) is a subgraph
of Ĝ. As a consequence, the statement of Theorem 4.2 can be strengthened: given any
inverse M-matrix Z 0 such that Z 0 ≥ R and Zii0 = Rii for all i, then ij ∈ Ĝ for all i 6= j
such that Zij0 = Rij .
Theorem 4.2 gives a partial justification for the heuristics used in applied statistics
to learn high-dimensional Gaussian graphical models by initializing search procedures
in MWSF(|R|), see for example [17]. We emphasize that the result is not true without
the MTP2 assumption.
An algebraic analogue of Theorem 4.2 holds for general inverse M-matrices. More
precisely, we have the following result:
Theorem 4.3. Assume that C is an inverse M-matrix with Cii = 1 and let K = C −1 .
Then if MWSF(C) consists of a single element, it holds that MWSF(C) ⊆ G(K) .
Proof. Let A denote the adjacency matrix of the complement of G(K) and let R =
C − εA. If ε is sufficiently small, we have MWSF(R) = MWSF(C) and Σ̂(R) = C
since C fulfills all the conditions in Proposition 2.2. Now we proceed as in the proof
of Theorem 4.2 and note that for uv 6∈ G(K), we have Cuv − Ruv = ε > 0 which
by complementary slackness implies Kuv = 0 and hence all edges in MWST(R) =
MWSF(C) must be in G(K).

Note that Theorem 4.2 and Theorem 4.3 represent a very strong form of absence of a
paradox of Yule–Simpson type for Gaussian MTP2 distributions. For MTP2 distributions, it clearly holds that a positive correlation cannot change sign when conditioning
with other variables. The results in these theorems are stronger and say that for Gaussian MTP2 distributions, correlations along edges of MWSF(C) necessarily correspond
to direct associations in the concentration graph and are not resulting from indirect
associations.
We have seen that the MWSF(R) provides a lower bound for the ML graph Ĝ =
G(K̂). Below we provide a simple procedure for identifying an upper bound for Ĝ.
This procedure relies on the estimation of the standard Gaussian graphical model over
the tree MWSF(R). The MLE under this assumption, denoted by Σ̃, can be computed
very efficiently and it satisfies
Y
Σ̃ij =
Ruv .
(u,v)∈ij

where ij denotes the unique path between i and j in MWSF(R); see, for example, [38,
Section 8.2].
Proposition 4.4. The pair ij forms an edge in the ML graph only if
Y
Rij ≥ Σ̃ij =
Ruv .
(u,v)∈ij

In particular Rij ≤ 0 implies that ij is not an edge of the ML graph.

16

MAXIMUM LIKELIHOOD IN TOTALLY POSITIVE GAUSSIAN MODELS

Proof. Because Σ̂ is an inverse M-matrix it is necessarily a path product matrix; see
[18, Section 1]. In particular, for all i, j
Y
Σ̂ij ≥
Σ̂uv .
(u,v)∈ij

By Corollary 2.4, the MLE under MTP2 is equal to the MLE under the Gaussian
graphical model over the ML graph. Therefore, if ij forms an edge in the ML graph,
then Σ̂ij = Rij . By Theorem 4.2, all edges of MWSF(R) lie in the ML graph and so
Σ̂uv = Ruv for all edges of ij. This proves that the inequality in the statement of the
proposition must hold.

Motivated by this result we define the excess correlation graph (EC graph) EC(R)
of R by the condition
Y
i ∼ j ⇐⇒ Rij ≥
Ruv .
(u,v)∈ij

Thus the EC graph has edges ij whenever the observed correlation between i and j is
in excess of or equal to what is explained by the spanning forest; by construction, we
have with probability one that
MWSF(R) ⊆ Ĝ = G(K̂) ⊆ EC(R).
The inclusions above are typically strict. For example, if R is an inverse M-matrix, then
EC(R) is the complete graph, whereas G(K̂) can be arbitrary; this follows from [12,
Proposition 6.3].
The following construction may be useful for studying MTP2 Gaussian distributions.
Let R be the sample correlation matrix. Define the matrix W by
Y
Wij := max
Ruv ,
P

(u,v)∈P

where, like in (12), the maximum is taken over all paths in G+ between i and j and is
set to zero if no such path exists. Transforming Dij = − log Rij gives a distance based
interpretation, in which Wij is related to the shortest distance between i and j in G+
with edge lengths given by Duv .
It is easy to show that Z ≥ W ≥ R and that W is always equal to the MLE Σ̂ in
the case when p ≤ 3. For general p we do not know conditions on R that assure that
W is an inverse M-matrix, or the MLE. We note that, quite surprisingly, W = Σ̂ for
the carcass data discussed in the introduction.
To complete the discussion of the ML graph, we provide a more detailed analysis in
the case when n = 2.
Proposition 4.5. Let R be a p × p sample correlation matrix of rank 2 such that
Rij < 1 for all i 6= j. Then the ML graph G(K̂) consists of at most two paths and thus
MWSF(R) = G(K̂).

MAXIMUM LIKELIHOOD IN TOTALLY POSITIVE GAUSSIAN MODELS

17

Proof. R can be decomposed into R = XX T , where X is a matrix of size p × 2
consisting of row vectors of length 1. The p row vectors of R are denoted by x1 , . . . , xp .
Without loss of generality we can assume that all vectors lie in the upper half-plane
and that x1 = (1, 0). In addition, without loss of generality we assume that the vectors
x1 , . . . , xp are arranged as shown in Figure 2. We denote by θij the angle between xi and
xj . Note that Rij = cos(θij ). Since cos(x) is decreasing on [0, π], the maximum weight
spanning tree is the path (1, 2, 3, . . . , p), assuming that all angles θi,i+1 ≤ π/2 (otherwise
the maximum weight spanning forest has two paths). By the previous theorem this
implies that these paths are contained in the ML graph. Buhl proved in [6] that the
MLE for the cycle with 2 observations exists if and only if the vectors x1 , . . . , xp do not
follow the same ordering as the cycle ordering. This implies that the MLE cannot be
a cycle. The result follows.

4
p

3
2
1

Figure 2. Configuration of vectors x1 , . . . , xp .
Note that since the ML graph is a forest in this case, the EC graph will consist of
complete graphs corresponding to the (at most two) connected components of the ML
graph.
5. Gaussian signed MTP2 distributions
In this section we discuss how our results can be generalized to so-called signed
MTP2 Gaussian distributions, where the distribution is MTP2 up to sign swapping.
Such distributions were discussed by Karlin and Rinott [21]. More precisely, a random
variable X has a signed MTP2 distribution if there exists a diagonal matrix D with
Dii = ±1 such that DX is MTP2 . In the context of Gaussian graphical models, signed
MTP2 distributions were studied in [26], where they are called non-frustrated.
The following characterization of signed MTP2 Gaussian distributions is a direct
consequence of [21, Theorem 3.1 and Remark 1.3].
Proposition 5.1. A Gaussian random variable X has a signed MTP2 distribution if
and only if |X| is MTP2 .
Our first result in this section is an alternative characterization of signed MTP2
Gaussian distributions via inequalities on the entries of the concentration matrix. The
proof is somewhat technical and can be found in the appendix.
Theorem 5.2. A Gaussian random variable X with concentration matrix K has a
signed MTP2 distribution if and only if for every cycle (i1 , . . . , ik , i1 ) in the graph G(K)
(13)

(−1)k Ki1 i2 Ki2 i3 · · · Kik i1 > 0.

18

MAXIMUM LIKELIHOOD IN TOTALLY POSITIVE GAUSSIAN MODELS

Signed MTP2 distributions are relevant for example because of their appearance
when studying tree models.
Proposition 5.3. Every Gaussian graphical model over a tree consists of signed MTP2
distributions. These MTP2 distributions are precisely the ones in which the covariance
matrix Σ has only nonnegative entries.
Proof. Let T be a tree and K = Σ−1 be a concentration matrix in the Gaussian
graphical model over T . Then G(K) is a subgraph of T and in particular it has no
cycles. Hence by Theorem 5.2 it is signed MTP2 . The second part of the statement
follows from [12, Corollary 7.3].

Because signed MTP2 distributions are closed under taking margins, Proposition 5.3
can be further generalized. The following theorem covers, in particular, Examples
4.1-4.5 in [21].
Theorem 5.4. Every distribution on a Gaussian tree model with hidden variables is
signed MTP2 .
Gaussian tree models with hidden variables have many applications, in partiular
related to modeling evolutionary processes; see, e.g., [7, 32]. As an important submodel
they contain the Brownian motion tree model [13]. Another example of a Gaussian tree
model is the factor analysis model with a single factor; it corresponds to a Gaussian
model on a star tree, whose inner node is hidden. The MTP2 distributions in this
model correspond to the distributions in a Spearman model [24, 34], where the hidden
factor is interpreted as intelligence.
Let R be a sample correlation matrix. Maximizing the likelihood over all signed
MTP2 Gaussian distributions corresponds to determining the sign matrix D, with Dii =
±1, that maximizes the likelihood for all 2p possible matrices DRD. A natural heuristic
∗
is to choose D∗ such that Dii∗ Djj
Rij ≥ 0 for all edges ij of MWSF(|R|). We provide
conditions under which this procedure indeed gives the MLE under signed MTP2 , and
we also provide examples showing that this is not true in general.
∗
Rij ≥ 0
First we describe how to obtain a sign swapping matrix D∗ such that Dii∗ Djj
for all edges ij of MWSF(|R|). Root MWSF(|R|) at node 1, that is, regard MWSF(|R|)
∗
as a directed tree with all edges directed away from 1. Set D11
= 1. Then proceed
∗
∗
recursively. For any edge i → j suppose that Dii is known and set Djj
:= sgn(Dii∗ Rij ).
Note that by construction

(14)

Dii∗ := sgn(R1i1 Ri1 i2 · · · Rik i ),

where 1 → i1 → i2 → · · · → ik → i is the unique path from 1 to i in MWSF(|R|). We
set Dii∗ = 0 if no such path exists. It is easy to check that the resulting D∗ satisfies
∗
Dii∗ Djj
Rij ≥ 0 for all edges ij of MWSF(|R|).
Proposition 5.5. Suppose that R is a sample correlation matrix such that
(15)

Rij Rik Rjk > 0 for all i, j, k ∈ V.

Then the MLE based on R over signed MTP2 Gaussian distributions is equal to the
MLE based on the sample correlation matrix D∗ RD∗ over MTP2 distributions.

MAXIMUM LIKELIHOOD IN TOTALLY POSITIVE GAUSSIAN MODELS

19

Proof. Let i1 − . . . − ik − i1 be an undirected cycle. Then (15) implies that
Ri1 i2 Ri2 i3 · · · Rik i1 > 0,

(16)

because the sign of this expression is equal to the sign of
(Ri1 i2 Ri2 i3 Ri1 i3 )(Ri1 i3 Ri3 i4 Ri1 i4 )(Ri1 i4 Ri4 i5 Ri1 i5 ) · · · (Ri1 ik Rik i1 Ri1 i1 ).
We now use this to show that D∗ RD∗ has only positive entries. Let i, j be any two
nodes and let 1 → i1 → · · · → ik → i and 1 → j1 → · · · → jl → j be the paths from 1
to i and j, respectively. By (14) we obtain
∗
Rij ) = sgn(R1i1 · · · Rik i Rij Rjjl · · · Rj1 1 ),
sgn(Dii∗ Djj

which is positive by (16).
So without loss of generality we can assume that all entries of R are nonnegative
and hence that D∗ is the identity matrix Ip . To complete the proof, we need to show
that the likelihood over MTP2 distributions given the sample correlation matrix DRD
is maximized by D = Ip . This is because (Dii Djj − 1) ≤ 0 and Rij Kij ≤ 0, and hence
X
`(K; R) − `(K; DRD) = tr(DRDK) − tr(RK) =
(Dii Djj − 1)Rij Kij ≥ 0,
i,j

which completes the proof.



By combining Proposition 5.5 with Proposition 4.5, we obtain the following result.
Corollary 5.6. If the sample correlation matrix R is based on n = 2 observations, then
the MLE over signed MTP2 Gaussian distributions given R is equal to the MLE over
MTP2 Gaussian distributions given the modified sample correlation matrix D∗ RD∗ and
consists of a single path.
Proposition 5.5 does not extend to an arbitrary sample correlation matrix. We give
a simple counterexample.
Example 5.7. Suppose that the sample correlation matrix is


1
0.3 0.11 0.3
 0.3
1
−0.1 −0.1
.
R=
0.11 −0.1
1
−0.1
0.3 −0.1 −0.1
1
Then MWSF(|R|) is given by the star graph with edges 1 − 2, 1 − 3, 1 − 4. Since R is
positive on these entries, D∗ = Ip . But one can check that the corresponding MLE has
a lower likelihood than the MLE after changing the sign of the third variable.
The intuition is the following. The log-likelihood based on R is up to an additive
constant given by
minimize − log det(Σ)
Σ

subject to Σ11 = Σ22 = Σ33 = Σ44 = 1,
Σ12 ≥ R12 , Σ13 ≥ R13 , Σ14 ≥ R14 ,
Σ23 ≥ 0, Σ24 ≥ 0, Σ34 ≥ 0,
Σ  0.

friendl

kind

approvn

sociabl
agreebl
coopera

talkatv
outgoin

hardwrk
discipl

organiz
persevr

laidbck
respnsi

relaxed

lax
easygon

withdrw

distant

quiet

shy

tense

anxious

worryin

givinup

lazy

carelss
disorgn

contrar

criticl
opposng

MAXIMUM LIKELIHOOD IN TOTALLY POSITIVE GAUSSIAN MODELS
harsh

20

1

harsh
criticl
opposng

0.8

contrar
carelss
disorgn

0.6

lazy
givinup
worryin
anxious

0.4

tense
shy
quiet

0.2

distant
withdrw
lax

0

easygon
relaxed
laidbck

−0.2

respnsi
organiz
persevr

−0.4

hardwrk
discipl
talkatv
outgoin

−0.6

sociabl
agreebl
coopera

−0.8

approvn
kind
friendl

−1

Figure 3. Correlation matrix of personality traits from the data set
described in [27].
By changing the sign of the third variable, we replace the constraint 1 − 3 by two
constraints 2 − 3 and 3 − 4. The resulting optimization problem is
minimize − log det(Σ)
Σ

subject to Σ11 = Σ22 = Σ33 = Σ44 = 1,
Σ12 ≥ R12 , Σ14 ≥ R14 , Σ23 ≥ −R23 , Σ34 ≥ −R34 ,
Σ13 ≥ 0, Σ24 ≥ 0,
Σ  0.
Note that R13 is only slightly larger than −R23 and −R24 . Hence, in essence we are
increasing the number of constraints by one, which explains the decrease of the loglikelihood value.
We end this paper by illustrating how our results can be applied to factor analysis
in psychometrics.
Example 5.8. Single factor models are routinely used to study the personalities in
psychometrics. Consider the following example from [27]3: 240 individuals were asked
to rate themselves on the scale 1-9 with respect to 32 different personality traits. The
3We

downloaded the
FactorAnalysis.html.

data

from

http://web.stanford.edu/class/psych253/tutorials/

worryin

1
givinup

anxious

persevr

tense

lazy

easygon

0.8
hardwrk

relaxed

discipl

laidbck

lax

harsh

carelss
disorgn

contrar

organiz

agreebl

harsh

coopera

criticl

approvn

opposng

distant
withdrw

0
approvn

shy

kind

relaxed

tense
easygon

anxious

sociabl
worryin

talkatv
outgoin

quiet

withdrw
shy

distant

friendl

kind

approvn

contrar
agreebl
coopera

criticl
opposng

harsh

organiz

carelss
disorgn

laidbck

1

0.8

0.6

0.4

0.2

0

friendl

talkatv

−0.2
distant

outgoin

withdrw

sociabl

shy

givinup

quiet
−0.4

persevr

talkatv

lazy

outgoin
sociabl
−0.6

discipl

worryin

lax

anxious

respnsi

lax
respnsi

0.2

contrar
agreebl
coopera

hardwrk

hardwrk
discipl

0.4

friendl

quiet

21

respnsi
0.6

criticl
opposng

kind

lazy

givinup
persevr

organiz

carelss
disorgn

lax
respnsi

hardwrk
discipl

lazy

sociabl
givinup
persevr

talkatv
outgoin

quiet

withdrw
shy

distant

friendl

kind

approvn

contrar
agreebl
coopera

criticl
opposng

harsh

laidbck

relaxed

tense
easygon

anxious

worryin

MAXIMUM LIKELIHOOD IN TOTALLY POSITIVE GAUSSIAN MODELS

tense

−0.8

carelss

easygon

disorgn

relaxed

organiz

laidbck

−1

−0.2

−0.4

−0.6

−0.8

−1

Figure 4. The correlation matrix of the data set on personality traits
after performing the sign switches as defined in (14) is shown on the
left. The correlation matrix resulting from switching the signs of the 16
(negative) traits that constitute the first block of variables in Figure 3 is
shown on the right.
resulting correlation matrix is shown in Figure 3. It appears to have a block structure
with predominantly positive entries in each diagonal block and negative entries in the
off-diagonal block. Also analyzing the respective variables, they seem to correspond to
positive and negative traits. It is therefore natural to assume that this data follows a
signed MTP2 distribution and analyze it under this constraint.
The correlation matrix resulting from the sign switching procedure described in (14)
is shown on the left in Figure 4, while the correlation matrix resulting from switching the
signs of the 16 (negative) traits that constitute the first block of variables in Figure 3 is
shown on the right in Figure 4. These plots suggest that the matrix on the right is closer
to being MTP2 . In fact, its log-likelihood (i.e., the value of n2 (log det K − tr(SK))) is
-2046.146, as compared to the log-likelihood value of -2071.717 resulting from the sign
switching procedure described in (14). For comparison, the value of the unconstrained
log-likelihood is -1725.075 and the value of the log-likelihood under MTP2 without sign
switching is -2356.639. The graphical models based on no sign switching and switching
the signs of the 16 negative traits are shown in Figure 5. The unconstrained loglikelihood gives a lower bound of 642.142 on the likelihood ratio statistic to test signed
MTP2 constraints, while the likelihood ratio statistic to test MTP2 constraints against
the saturated model is equal to 1263.128.
Appendix A. Proof of Theorem 5.2
We start by proving the following auxiliary result.
Proposition A.1. A (non-degenerate) Gaussian random variable X has a signed
MTP2 distribution if and only if its concentration matrix lies in the closure of the

22

MAXIMUM LIKELIHOOD IN TOTALLY POSITIVE GAUSSIAN MODELS
13

15

15

17
4

4

12

17

8

22

25

2

23

20

24

6

13
26

18
28

25
31

28

22

8
2

23

9

24
21

1

27

32

7

3

26

6

5

11

30

20

12

29

10

31

30
32

14

1

14
18

16

29

16

10

21
19
7

19

9
27

5

Figure 5. The graphical models resulting from estimation under MTP2
based on the correlation matrix shown in Figure 3 and the correlation
matrix shown in Figure 4 (right). The thin gray edges correspond to the
edges of the EC graph that are not part of the ML graph. The blue
edges represent edges of the ML graph that are not part of the minimum
weight spanning tree. The latter is represented by thick red edges.
set satisfying
(17)

Kij Kik Kjk < 0

for all 1 ≤ i < j < k ≤ p.

Proof. If K̃ := DKD is an M-matrix then
2
2
K̃ij K̃ik K̃jk = Dii2 Djj
Dkk
Kij Kik Kjk = Kij Kik Kjk

and so Kij Kik Kjk ≤ 0 for all 1 ≤ i < j < k ≤ p . If K̃ has no zeros, then K lies in the set
defined by (17) and so also in its closure. If K̃ has zeros, there exists a sequence K̃ (n) of
matrices with strictly negative off-diagonal elements such that K̃ (n) → K̃ elementwise.
In that case also DK̃ (n) D → DK̃D = K and so K lies in the closure of the set defined
by (17). This proves the ’only if’ direction.
Now suppose that K lies in the closure of the set given by (17). First, consider the
case, when K satisfies (17) and thus, in particular, Kij 6= 0 for all i < j. Then
K̃ij := (sgnK1i sgnK1j )Kij = −(sgnKij )Kij = −|Kij | < 0,
where the second equation follows from the fact that K1i K1j Kij < 0. In particular,
K̃ = [K̃ij ] is an M-matrix. If some Kij are zero, then we take a sequence K (n) of
matrices satisfying (17) such that K (n) → K. For each K (n) we determine D(n) such
that D(n) K (n) D(n) is an M-matrix. By continuity of matrix multiplication and the fact

11

3

MAXIMUM LIKELIHOOD IN TOTALLY POSITIVE GAUSSIAN MODELS

23

that the set of M-matrices is relatively closed, there must exist D such that DKD is
an M-matrix and D = lim D(n) .

Remark A.2. It may look surprising at first but the closure of the set given by (17) is
not the set with strict inequalities replaced by ≤. For example, consider a 4 × 4-matrix
K such that K13 = K24 = 0, K12 > 0, and K23 , K34 , K14 < 0. Then for all triples it
(n)
holds that Kij Kik Kjk ≤ 0 (in fact “= 0”), but there is no sequence K13 → 0 such that
(n)
(n)
K12 K13 K23 < 0 and K13 K14 K34 < 0.
Proof of Theorem 5.2. To prove necessity, note that every point that lies in the closure
of (17) satisfies also the non-strict version of the inequalities in (17). The inequalities
Kij Kik Kjk ≤ 0 for triples {i1 , i2 , i3 }, {i1 , i3 , i4 }, {i1 , i4 , i5 }, . . . , {i1 , ik−1 , ik } multiplied
together yield (−1)k Ki1 i2 Ki2 i3 · · · Kik i1 ≥ 0. Because all elements in this product are
nonzero, (13) must hold.
To prove sufficiency, by Proposition A.1 it suffices to show that every point K satisfying (13) lies in the closure of the set described by (17). If there is a vertex, say
vertex 1, in G(K) that is linked by an edge to any other vertex, we can proceed as in
the proof of Proposition A.1. Let D be a diagonal matrix such that Dii = −sgnK1i .
Let K̃ = DKD. If Kij = 0, then K̃ij = 0. If Kij 6= 0, then
K̃ij = Dii Djj Kij = sgn(K1i )sgn(K1j )Kij = −|Kij |,
where the last equality follows from the fact that necessarily K1i K1j Kij < 0. Hence
K̃ is an M-matrix and thus K is the concentration matrix of a signed MTP2 Gaussian
distribution. As a consequence, it lies in the closure of (17).
Now suppose that G(K) does not contain a vertex linked by an edge to any other
vertex. If G(K) is not connected we can use an inductive argument because K must
be block-diagonal. Hence we can assume that G(K) is connected. To show that K lies
(n)
in the closure of (17) we construct a sequence K (n) → K as follows. Let Kij = Kij
(n)
whenever i, j 6= 1. If K1j 6= 0 then let K1j = K1j for all n ≥ 1. If K1j = 0 then
consider any path (1, i1 , . . . , ik , j) from 1 to j in G(K). The sign s ∈ {−1, +1} of the
expression (−1)k K1i1 Ki1 i2 · · · Kik j does not depend on the choice of the path. Indeed,
if (1, j1 , . . . , jl , j) is any other path, then the walk (1, i1 , . . . , ik , j, jl , . . . , j1 , 1) can be
decomposed into cycles in G(K). For each of these cycles we have inequality (13).
Multiplying all these inequalities together yields
(−1)k K1i1 Ki1 i2 · · · Kik j · (−1)l Kjjl · · · Kj2 j1 ≥ 0.
Hence the signs of the two paths are either both negative or both positive. Now for all
(n)
(n)
j such that K1j = 0 set K1j = − ns . By construction, K (n) → K and K1j 6= 0 for all
j ∈ V and n ≥ 1. Then as in the previous paragraph we construct a diagonal matrix
(n)
(n)
(n)
D(n) such that Dii = −sgn(K1i ) and let K̃ (n) = D(n) K (n) D(n) . If Kij = 0, then
(n)
(n)
K̃ij = 0. If Kij 6= 0, then
(n)

(n)

(n)

(n)

(n)

K̃ij = sgn(K1i )sgn(K1j )Kij = −|Kij |.
(n)

(n)

(n)

The proof of the last equality is based on the fact that K1i K1j Kij < 0, which we
(n) (n) (n)
now show to conclude the proof. By construction, the sign of K1i K1j Kij is equal to

24

MAXIMUM LIKELIHOOD IN TOTALLY POSITIVE GAUSSIAN MODELS

the sign of
(18)

(−1)k+l K1i1 · · · Kik i Kij Kjjl · · · Ki1 1 ,

where (1, i1 , . . . , ik , i) and (1, j1 , . . . , jl , j) are any two paths in GK ; take (1, i) if K1i 6= 0,
and (1, j) if K1j 6= 0. By splitting the walk (1, i2 , . . . , ik , i, j, jl , . . . , i1 , 1) into cycles, we
can use (13) to show that (18) must be negative.

References
[1] A. Anandkumar, V. Y. Tan, F. Huang, and A. S. Willsky, High-dimensional Gaussian graphical model selection: Walk summability and local separation criterion, The Journal of
Machine Learning Research, 13 (2012), pp. 2293–2337.
[2] F. Bartolucci and J. Besag, A recursive algorithm for Markov random fields, Biometrika, 89
(2002), pp. 724–730.
[3] F. Bartolucci and A. Forcina, A likelihood ratio test for MTP2 within binary variables, Ann.
Statist., 28 (2000), pp. 1206–1218.
[4] B. Bhattacharya, Covariance selection and multivariate dependence, Journal of Multivariate
Analysis, 106 (2012), pp. 212–228.
[5] E. Bølviken, Probability inequalities for the multivariate normal with non-negative partial correlations, Scand. J. Statist., 9 (1982), pp. 49–58.
[6] S. L. Buhl, On the existence of maximum likelihood estimators for graphical Gaussian models,
Scand. J. Statist., 20 (1993), pp. 263–270.
[7] M. J. Choi, V. Y. F. Tan, A. Anandkumar, and A. S. Willsky, Learning latent tree
graphical models, J. Mach. Learn. Res., 12 (2011), pp. 1771–1812.
[8] A. Colangelo, M. Scarsini, and M. Shaked, Some notions of multivariate positive dependence, Insurance: Mathematics and Economics, 37 (2005), pp. 13–26.
[9] C. Dellacherie, S. Martinez, and J. San Martin, Inverse M-matrices and ultrametric
matrices, vol. 2118, Springer, 2014.
[10] A. P. Dempster, Covariance selection, Biometrics, (1972), pp. 157–175.
[11] J. Djolonga and A. Krause, Scalable variational inference in log-supermodular models, In
International Conference on Machine Learning (ICML), (2015).
[12] S. Fallat, S. L. Lauritzen, K. Sadeghi, C. Uhler, N. Wermuth, and P. Zwiernik,
Total positivity in Markov structures, to appear in the Annals of Statistics, (2016).
[13] J. Felsenstein, Maximum-likelihood estimation of evolutionary trees from continuous characters., American Journal of Human Genetics, 25 (1973), pp. 471–492.
[14] C. M. Fortuin, P. W. Kasteleyn, and J. Ginibre, Correlation inequalities on some partially
ordered sets, Comm. Math. Phys., 22 (1971), pp. 89–103.
[15] J. C. Gower and G. Ross, Minimum spanning trees and single linkage cluster analysis, Applied
statistics, (1969), pp. 54–64.
[16] M. Grant and S. Boyd, CVX: Matlab software for disciplined convex programming, version
2.1. http://cvxr.com/cvx, Mar. 2014.
[17] S. Højsgaard, D. Edwards, and S. Lauritzen, Graphical Models with R, Springer, New
York, 2012.
[18] C. R. Johnson and R. L. Smith, Path product matrices, Linear and Multilinear Algebra, 46
(1999), pp. 177–191.
[19] C. R. Johnson and R. L. Smith, Inverse M-matrices, II, Linear Algebra and its Applications,
435 (2011), pp. 953–983.
[20] S. Karlin and Y. Rinott, Classes of orderings of measures and related correlation inequalities.
I. Multivariate totally positive distributions, J. Multiv. Anal., 10 (1980), pp. 467–498.
, Total positivity properties of absolute value multinormal variables with applications to con[21]
fidence interval estimates and related probabilistic inequalities, Ann. Statist., 9 (1981), pp. 1035–
1049.

MAXIMUM LIKELIHOOD IN TOTALLY POSITIVE GAUSSIAN MODELS

[22]
[23]
[24]
[25]

[26]
[27]
[28]
[29]
[30]
[31]
[32]
[33]

[34]
[35]
[36]
[37]

[38]

25

, M-matrices as covariance matrices of multinormal distributions, Linear Algebra and its
Applications, 52 (1983), pp. 419 – 438.
S. L. Lauritzen, Graphical Models, Clarendon Press, Oxford, United Kingdom, 1996.
W. Ledermann, I.—On a problem concerning matrices with variable diagonal elements., Proceedings of the Royal Society of Edinburgh, 60 (1940), pp. 1–17.
Z. Q. Luo and P. Tseng, On the convergence of the coordinate descent method for convex
differentiable minimization, Journal of Optimization Theory and Applications, 72 (1992), pp. 7–
35.
D. M. Malioutov, J. K. Johnson, and A. S. Willsky, Walk-sums and belief propagation in
Gaussian graphical models, The Journal of Machine Learning Research, 7 (2006), pp. 2031–2064.
B. F. Malle and L. M. Horowitz, The puzzle of negative self-views: An exploration using the
schema concept., Journal of Personality and Social Psychology, 68 (1995), p. 470.
C. M. Newman, A general central limit theorem for FKG systems, Comm. Math. Phys., 91
(1983), pp. 75–80.
A. Ostrowski, Über die Determinanten mit überwiegender Hauptdiagonale, Commentarii Mathematici Helvetici, 10 (1937), pp. 69–96.
J. G. Propp and D. B. Wilson, Exact sampling with coupled Markov chains and applications
to statistical mechanics, Random Structures and Algorithms, 9 (1996), pp. 223–252.
C. Semple and M. A. Steel, Phylogenetics, vol. 24, Oxford University Press, 2003.
N. Shiers, P. Zwiernik, J. Aston, and J. Q. Smith, The correlation space of Gaussian latent
tree models and model selection without fitting, Biometrika, 103(3) (2016), pp. 531–545.
M. Slawski and M. Hein, Estimation of positive definite M-matrices and structure learning
for attractive Gaussian Markov random fields, Linear Algebra and its Applications, 473 (2015),
pp. 145–179.
C. Spearman, The abilities of man, Science, 68 (1928), p. 38.
T. P. Speed and H. T. Kiiveri, Gaussian Markov distributions over finite graphs, Ann. Statist.,
14 (1986), pp. 138–150.
C. Uhler, Geometry of maximum likelihood estimation in Gaussian graphical models, Ann.
Statist., 40 (2012), pp. 238–261.
N. Wermuth and E. Scheidt, Algorithm AS 105: Fitting a covariance selection model to a
matrix, Journal of the Royal Statistical Society. Series C (Applied Statistics), 26 (1977), pp. pp.
88–92.
P. Zwiernik, Semialgebraic Statistics and Latent Tree Models, no. 146 in Monographs on Statistics and Applied Probability, Chapman & Hall, 2015.

Department of Mathematical Sciences, University of Copenhagen, Copenhagen, Denmark
E-mail address: lauritzen@math.ku.dk
Laboratory for Information and Decision Systems, and Institute for Data, Systems
and Society, Massachusetts Institute of Technology, Cambridge, MA, USA
E-mail address: cuhler@mit.edu
Department of Economics and Business, Universitat Pompeu Fabra, Barcelona, Spain
E-mail address: piotr.zwiernik@upf.edu

