Maximum likelihood estimation in Gaussian models under total positivity

TL;DR

This paper studies maximum likelihood estimation for Gaussian models with total positivity constraints, proving existence with minimal data, and introduces algorithms for efficient estimation and analysis of the resulting sparse inverse covariance structures.

Contribution

It provides a simple proof of MLE existence for MTP2 Gaussian models with minimal observations and develops globally convergent algorithms for estimation under these constraints.

Findings

MLE exists with at least 2 observations regardless of dimension

The ML graph can be bounded using maximum weight spanning forest

Algorithms similar to iterative proportional scaling are proposed

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 at least 2 observations, irrespective of the underlying dimension. Slawski and Hein, 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 we can find an upper bound for the ML graph by adding edges corresponding to correlations in excess of those explained by the maximum weight spanning forest of the correlation matrix. Moreover, 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.

Tables

Table 1. Table 1: Vertex labeling for Figure 4 .

1	2	3	4	5	6	7	8
distant	talkatv	carelss	hardwrk	anxious	agreebl	tense	kind
9	10	11	12	13	14	15	16
opposng	relaxed	disorgn	outgoin	approvn	shy	discipl	harsh
17	18	19	20	21	22	23	24
persevr	friendl	worryin	respnsi	contrar	sociabl	lazy	coopera
25	26	27	28	29	30	31	32
quiet	organiz	criticl	lax	laidbck	withdrw	givinup	easygon

Equations

$$f(x)f(y)\le f(x\wedge y)f(x\vee y).$$

$$\ell (K;S)=\log \det K-\mathrm{tr}(SK).$$

$$\small R=\blockarray{ccccccc}\textrm{Fat11}&\textrm{Meat11}\textrm{Fat12}\textrm{Meat12}\textrm{Fat13}\textrm{Meat13}\\ \block{(rrrrrr)l}1.000.040.840.080.82-0.03\textrm{Fat11}\\ 0.041.000.040.870.130.86\textrm{Meat11}\\ 0.840.041.000.010.83-0.03\textrm{Fat12}\\ 0.080.870.011.000.110.90\textrm{Meat12}\\ 0.820.130.830.111.000.02\textrm{Fat13}\\ -0.030.86-0.030.900.021.00\textrm{Meat13}\\$$

$$\small\tilde{K}=\blockarray{ccccccc}\textrm{Fat11}&\textrm{Meat11}\textrm{Fat12}\textrm{Meat12}\textrm{Fat13}\textrm{Meat13}\\ \block{(rrrrrr)l}1.00{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0.16}-0.52-0.31-0.40{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0.19}\textrm{Fat11}\\ {\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0.16}1.00-0.05-0.42-0.17-0.37\textrm{Meat11}\\ -0.52-0.051.00{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0.25}-0.45-0.17\textrm{Fat12}\\ -0.31-0.42{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0.25}1.00-0.02-0.61\textrm{Meat12}\\ -0.40-0.17-0.45-0.021.00{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0.10}\textrm{Fat13}\\ {\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0.19}-0.37-0.17-0.61{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0.10}1.00\textrm{Meat13}\\$$

$$\small\hat{R}=\blockarray{ccccccc}\textrm{Fat11}&\textrm{Meat11}\textrm{Fat12}\textrm{Meat12}\textrm{Fat13}\textrm{Meat13}\\ \block{(rrrrrr)l}1.00{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.10}0.84{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.09}0.82{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.09}\textrm{Fat11}\\ {\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.10}1.00{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.11}0.870.130.86\textrm{Meat11}\\ 0.84{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.11}1.00{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.09}0.83{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.09}\textrm{Fat12}\\ {\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.09}0.87{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.09}1.00{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.11}0.90\textrm{Meat12}\\ 0.820.130.83{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.11}1.00{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.11}\textrm{Fat13}\\ {\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.09}0.86{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.09}0.90{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.11}1.00\textrm{Meat13}\\$$

$$x\wedge y=(\min (x_v,y_v),v\in V),x\vee y=(\max (x_i,y_i),i\in V).$$

$$f(x)f(y)\le f(x\wedge y)f(x\vee y)\text{mbox}forallx,y\in \mathcal{X}.$$

$${\mathcal{H}}_{ij}^p=\{X\in {\mathbb{S}}^p\mid X_{ij}\le 0\}$$

$$\underset{K}{\mathrm{maximize}}$$

$$K\in {\mathcal{M}}^p$$

$${\mathcal{N}}^p=\{X\in {\mathbb{S}}^p\mid \exists Y\in {\mathbb{S}}_{\succ 0}^p\text{ with }X\le Y\text{ and }\text{diag}(X)=\text{diag}(Y)\}.$$

$$\overline{\mathcal{N}}\,=\,\bigl{\{}\,S\in\mathbb{S}\,\mid\,\langle S,K\rangle\geq 0\,\,\,\hbox{for all}\,\,\,K\in\mathcal{M}\bigr{\}}.$$

$$({\mathcal{C}}_1\cap {\mathcal{C}}_2{)}^{\vee }={\mathcal{C}}_1^{\vee }+{\mathcal{C}}_2^{\vee };$$

$${\mathbb{S}}_{\succ 0}^{\vee }={\mathbb{S}}_{⪰0}\text{ and }{\mathcal{H}}_{ij}^{\vee }={\mathcal{H}}_{ij}.$$

$$({\hat{\Sigma }}^{-1}{)}_{ij}$$

$${\hat{\Sigma }}_{ii}-S_{ii}$$

$$({\hat{\Sigma }}_{ij}-S_{ij})$$

$$({\hat{\Sigma }}_{ij}-S_{ij})({\hat{\Sigma }}^{-1}{)}_{ij}$$

$$\underset{\Sigma ⪰0}{\mathrm{minimize}}$$

$$\displaystyle\Sigma_{ii}=S_{ii},\;\mbox{ for all $i\in V$},$$

$$\displaystyle\Sigma_{ij}\geq S_{ij},\;\mbox{ for all $i\neq j$}.$$

$${\bar{\Sigma }}_{ii}-S_{ii}$$

$${\bar{\Sigma }}_{ij}-S_{ij}$$

$$({\bar{\Sigma }}^{-1}{)}_{ij}$$

$${\hat{\Sigma }}_{ij}^S=\sqrt{S_{ii}{S}_{jj}}{\hat{\Sigma }}_{ij}^R\text{mbox}foralli,j\in V.$$

$$\log \det K-\mathrm{tr}(SK)=\log \det K-\mathrm{tr}(RDKD).$$

$$\log \det K-\text{trace}(KS)=\log \det K^{\prime }-\text{trace}(K^{\prime }{S}_{AA}).$$

$$\underset{K^{\prime }⪰0}{\mathrm{maximize}}$$

$$K_{12}^{\prime }+L_{12}\le 0.$$

$$L_{12}\le \frac{S_{uv}}{S_{uu}{S}_{vv}-S_{uv}^2}.$$

$$K_{11}^{\prime }=\frac{1+\sqrt{1+4S_{uu}{S}_{vv}{L}_{12}^2}}{2S_{uu}},K_{22}^{\prime }=\frac{1+\sqrt{1+4S_{uu}{S}_{vv}{L}_{12}^2}}{2S_{vv}},K_{12}^{\prime }=-L_{12}.$$

$${\Sigma }_{uv}=\max \{S_{uv},L_{12}\}.$$

$$\underset{K⪰0}{\mathrm{maximize}}$$

$$K_{uv}\le 0,$$

$$K_{ij}=K_{ij}^1\text{ for all }ij\in (V\times V)\setminus \{uu,vv,uv\},$$

$$\underset{\Sigma ⪰0}{\mathrm{maximize}}$$

Full text

60E15,

62H99,

15B48,

${\rm MTP}_{2}$ distribution; attractive Gaussian Markov random field (GMRF); non-frustrated GRMF; Gaussian graphical model; inverse M-matrix; ultrametric.,

, , and

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, 23]. This property (also called the ${\rm MTP}_{2}$ property) appeared in the study of stochastic orderings, asymptotic statistics, and in statistical physics [15, 31]. Families of distributions with this property lead to many computational advantages [2, 11, 33]. In a recent paper [13], the ${\rm MTP}_{2}$ property was studied in the context of graphical models and conditional independence in general. It was shown that ${\rm MTP}_{2}$ 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 $\mu$ is multivariate totally positive of order 2 (${\rm MTP}_{2}$) if the density satisfies

[TABLE]

A multivariate Gaussian distribution with mean $\mu$ and a positive definite covariance matrix $\Sigma$ is ${\rm MTP}_{2}$ if and only if the concentration matrix $K:=\Sigma^{-1}$ is a symmetric M-matrix, that is, $K_{ij}\leq 0$ for all $i\neq j$ or, equivalently, if all partial correlations are nonnegative. Such distributions were considered by Bølviken [5] and Karlin and Rinott [25]. Moreover, Gaussian graphical models, or Gaussian Markov random fields, were studied in the context of totally positive distributions in [29]. ${\rm MTP}_{2}$ 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 ${\rm MTP}_{2}$ 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 [37] and this paper is closely related to their approach; see also [4] and [12].

Throughout this paper, we assume that we are given $n$ i.i.d. samples from $\mathcal{N}(\mu,\Sigma)$, where $\Sigma$ is an unknown positive definite matrix of size $p\times p$. Without loss of generality, we assume that $\mu=0$ and we focus on the estimation of $\Sigma$. We denote the sample covariance matrix based on $n$ samples by $S$. Then the log-likelihood function is, up to additive and multiplicative constants, given by

[TABLE]

We denote the cone of real symmetric matrices of size $p\times p$ by $\mathbb{S}^{p}$, its positive definite elements by $\mathbb{S}^{p}_{\succ 0}$, and its positive semidefinite elements by $\mathbb{S}^{p}_{\succeq 0}$. Note that $\ell(K;S)$ is a strictly concave function of $K\in\mathbb{S}^{p}_{\succeq 0}$ . Since M-matrices form a convex subset of $\mathbb{S}^{p}_{\succ 0}$, the optimization problem for computing the maximum likelihood estimator (MLE) for ${\rm MTP}_{2}$ Gaussian models is a convex optimization problem. Slawski and Hein [37] showed that the MLE exists, i.e., the global maximum of this optimization problem is attained, when $n\geq 2$. This yields a drastic reduction from $n\geq p$ without the ${\rm MTP}_{2}$ constraint. In addition, they provided empirical evidence showing that the ${\rm MTP}_{2}$ 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 $\hat{K}$ under the ${\rm MTP}_{2}$ constraint. For a $p\times p$ matrix $K$ we let $G(K)$ denote the undirected graph on $p$ nodes with an edge $ij$ if and only if $K_{ij}\neq 0$. In Proposition 4.3 we obtain a simple upper bound for the ML graph $G(\hat{K})$ by adding edges to the smallest maximum weight spanning forest (MWSF) corresponding to empirical correlations in excess of those provided by the MWSF. We illustrate the problem in the following example.

Example 1.1**.**

We consider the carcass data that are discussed in [19] 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 ignore the lean meat percentage, since, by definition, this variable should be negatively correlated with fat and positively correlated with meat so the joint distribution is unlikely to be ${\rm MTP}_{2}$. The sample correlation matrix $R$ for these data is

[TABLE]

and its inverse, scaled to have diagonal elements equal to one, $\tilde{K}$, is

[TABLE]

Note that the off-diagonal entries of $\tilde{K}$ are the negative empirical partial correlations. This sample distribution is not ${\rm MTP}_{2}$; the positive entries in $\tilde{K}$ are highlighted in red. The MLE under ${\rm MTP}_{2}$ can be computed for example using cvx [17] in matlab or using one of the simple coordinate descent algorithms discussed in Section 2. In this particular example the MLE can also be obtained through the explicit formula (14) in Section 4. The MLE of the correlation matrix, rounded to 2 decimals, is

[TABLE]

The entries of $\hat{R}$ that changed compared to the sample correlation matrix $R$ are highlighted in blue111We note that $\hat{\Sigma}_{45}>S_{45}$; the entries appear equal only because of the 2-digit rounding.. The sparsity pattern of $\hat{K}=\hat{\Sigma}^{-1}$ is captured by the ML graph $G(\hat{K})$ shown in Figure 1.

Note that all edges corresponding to blue entries in $\hat{R}$ 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 $R$. In this example, the spanning forest is a chain represented by the thick red edges in Figure 1. By Corollary 4.7 these edges form a spanning tree of the ML graph $G(\hat{K})$.

Interestingly, applying various methods for model selection such as stepwise AIC, BIC, or graphical lasso all yield similar graphs, possibly indicating that the ${\rm MTP}_{2}$ assumption is quite reasonable.

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 ${\rm MTP}_{2}$ Gaussian distributions. This embeds the results by Slawski and Hein [37] 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 $\Sigma$. In Section 3, we show how the problem of ML estimation for ${\rm MTP}_{2}$ 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\geq 2$, a result that was first proven in [37]. Our proof is by constructing a primal and dual feasible point of the convex ML estimation problem for ${\rm MTP}_{2}$ Gaussian models. In Section 4 we investigate the structure of the ML graph $G(\hat{K})$ and give a simple upper bound for it. Finally, in Section 5 we discuss how our results can be generalized to so-called signed ${\rm MTP}_{2}$ Gaussian distributions, where the distribution is ${\rm MTP}_{2}$ up to sign changes or, equivalently, $|X|$ is ${\rm MTP}_{2}$. Such distributions were introduced by Karlin and Rinott in [24]. We conclude the paper with a brief discussion of various open problems.

2 Duality theory for ML estimation under ${\rm MTP}_{2}$

We start this section by formally introducing absolutely continuous ${\rm MTP}_{2}$ distributions and then discuss the duality theory for Gaussian ${\rm MTP}_{2}$ distributions. Let $V:=\{1,2,\dots,p\}$ be a finite set and let $X=(X_{i},i\in V)$ be a random vector with density $f$ w.r.t. Lebesgue measure on the product space $\mathcal{X}=\prod_{i\in V}\mathcal{X}_{i}$, where $\mathcal{X}_{i}\subseteq\mathbb{R}$ is the state space of $X_{i}$. We define the coordinate-wise minimum and maximum as

[TABLE]

Then we say that $X$ or the distribution of $X$ is multivariate totally positive of order two (${\rm MTP}_{2}$) if its density function $f$ on $\mathcal{X}$ satisfies

[TABLE]

In this paper, we concentrate on the Gaussian setting. It is easy to show that a Gaussian distribution with mean $\mu$ and covariance matrix $\Sigma$ is ${\rm MTP}_{2}$ if and only if $K=\Sigma^{-1}$ is a symmetric M-matrix, i.e. $K$ is positive definite and

• (i)

$K_{ii}>0$ for all $i\in V$,

• (ii)

$K_{ij}\leq 0$ for all $i,j\in V$ with $i\neq j$.

Properties of M-matrices were studied by Ostrowski [32] who chose the name to honor H. Minkowski. The connection to multivariate Gaussian distributions was established by Bølviken [5] and Karlin and Rinott [25].

We denote the set of all symmetric M-matrices of size $p\times p$ by $\mathcal{M}^{p}$. Note that $\mathcal{M}^{p}$ is a convex cone. In fact, it is obtained by intersecting the positive definite cone $\mathbb{S}^{p}_{\succ 0}$ with all the coordinate half-spaces

[TABLE]

with $i\neq j$. For a convex cone $\mathcal{C}$ we denote its closure by $\overline{\mathcal{C}}$. Then $\overline{\mathcal{M}^{p}}$ is given by $\mathbb{S}^{p}_{\succeq 0}\cap_{i<j}\mathcal{H}^{p}_{ij}$ and the ML estimation problem for Gaussian ${\rm MTP}_{2}$ models can be formulated as the following optimization problem:

[TABLE]

This is a convex optimization problem, since the objective function is concave on $\mathbb{S}^{p}_{\succeq 0}$.

Next, we introduce a second convex cone $\mathcal{N}^{p}$ that plays an important role for ML estimation in Gaussian ${\rm MTP}_{2}$ models. To formally define this cone, we introduce two partial orders on matrices. Let $A,B$ be two $p\times p$ matrices. Then $A\geq B$ means that $A_{ij}\geq B_{ij}$ for all $(i,j)\in V\times V$, and $A\succeq B$ means that $A-B\in\mathbb{S}^{p}_{\succeq 0}$. Then the cone $\mathcal{N}^{p}$ is defined as the negative closure of $\mathbb{S}^{p}_{\succ 0}$, i.e.

[TABLE]

To simplify notation, we will suppress the dependence on $p$ and write $\mathbb{S}$, $\mathbb{S}_{\succeq 0}$, $\mathbb{S}_{\succ 0}$, $\mathcal{M}$ and $\mathcal{N}$, when the dimension is clear. In the following result, we show that the cones $\mathcal{N}$ and $\mathcal{M}$ are dual to each other.

Lemma 2.1**.**

The closure of $\mathcal{N}$ is the dual to the cone of M-matrices $\mathcal{M}$, i.e.

[TABLE]

Proof.

We denote the dual of a convex cone $\mathcal{C}$ by $\mathcal{C}^{\vee}$. Let $\mathcal{C}_{1}$, $\mathcal{C}_{2}$ be two convex cones. Then it is an easy exercise to verify that

[TABLE]

here $+$ denotes the Minkowski sum. Note that

[TABLE]

This completes the proof, since $\mathcal{M}=\mathbb{S}_{\succ 0}\cap_{i<j}\mathcal{H}_{ij}$ and (5) can be applied inductively to any finite collection of convex cones. ∎

Using the cones $\mathcal{M}$ and $\mathcal{N}$ we now determine conditions for existence of the MLE in Gaussian ${\rm MTP}_{2}$ 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 ${\rm MTP}_{2}$ model. Then the MLE $\hat{\Sigma}$ (and $\hat{K}$) exists for a given sample covariance matrix $S$ on $V$ if and only if $S\in\mathcal{N}$. It is then equal to the unique element $\hat{\Sigma}\succ 0$ that satisfies the following system of equations and inequalities

[TABLE]

Proof.

It is straight-forward to compute the dual optimization problem and the KKT conditions. In particular, in [37] it was shown that the dual optimization problem to (3) is given by

[TABLE]

Note that the identity matrix is a strictly feasible point for (3). Hence, the MLE does not exist if and only if the likelihood is unbounded. Since by Slater’s constraint qualification strong duality holds for the optimization problems (3) and (10), the MLE does not exist if and only if $S\notin\mathcal{N}$. ∎

We note that the conditions in Proposition 2.2 were also derived in [37], save for the explicit identification of the dual cone $\mathcal{N}$.

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 ${\rm MTP}_{2}$. In this case, let $G=(V,E)$ be an undirected graph. Then the primal problem has additional equality constraints, namely $K_{ij}=0$ for all $ij\notin E$, and hence the inequality constraints in the dual problem are restricted to the entries in $E$, i.e., $\Sigma_{ij}\geq S_{ij}$ for all $ij\in E$. Note that if the MLE of $\Sigma$ based on $S$ exists in the Gaussian graphical model over $G$, it also exists in the Gaussian graphical model over $G$ under ${\rm MTP}_{2}$, since without the ${\rm MTP}_{2}$ constraint the MLE needs to satisfy $\hat{\Sigma}_{ij}=S_{ij}$ for all $ij\in E$. ∎

We define the maximum likelihood graph (ML graph) $\hat{G}$ to be the graph determined by $\hat{K}$, i.e. $\hat{G}=G(\hat{K})$, where $\hat{K}=\hat{\Sigma}^{-1}$ is the MLE of $K$ under ${\rm MTP}_{2}$. We then have the following important corollary of Proposition 2.2.

Corollary 2.4**.**

Consider the Gaussian graphical model determined by $K_{ij}=0$ for $ij\not\in E(\hat{G})$, where $\hat{G}$ is the ML graph under ${\rm MTP}_{2}$. Let $\bar{\Sigma}$ be the MLE of $\Sigma$ under that Gaussian graphical model (without the ${\rm MTP}_{2}$ constraint). Then $\hat{\Sigma}=\bar{\Sigma}$.

Proof.

The MLE of $\Sigma$ under the Gaussian graphical model with graph $\hat{G}$ is the unique element $\bar{\Sigma}\succ 0$ satisfying the following system of equations:

[TABLE]

Proposition 2.2 says that also $\hat{\Sigma}$ satisfies these equations and hence we must have $\bar{\Sigma}=\hat{\Sigma}$. ∎

Note that this corollary highlights the role of the complementary slackness condition (9) in inducing sparsity of the ${\rm MTP}_{2}$ solution.

We emphasize that the MLE under ${\rm MTP}_{2}$ is equivariant w.r.t. changes of scale so that without loss of generality we can assume that the sample covariance is normalized, i.e. $S_{ii}=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 $\hat{\Sigma}^{S}$ and $\hat{\Sigma}^{R}$ the MLE in Proposition 2.2 based on $S$ and $R$, respectively. Then

[TABLE]

Proof.

Denote by $D$ a diagonal matrix such that $D_{ii}=\sqrt{S_{ii}}$ and $S=DRD$. The likelihood function based on $S$ is

[TABLE]

If $K^{\prime}=DKD$, this can be rewritten as $\log\det K^{\prime}-{\rm tr}(RK^{\prime})-\sum_{i}\log S_{ii}$. Therefore, if $\hat{K}^{R}$ is the maximizer of $\log\det K-{\rm tr}(RK)$ under the ${\rm MTP}_{2}$ constraints, then $D^{-1}\hat{K}^{R}D^{-1}$ is also an M-matrix and the maximizer of $\log\det K-{\rm tr}(SK)$. ∎

We end this section by providing simple coordinate descent algorithms for ML estimation under ${\rm MTP}_{2}$. 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\times 2$ which have a simple and explicit solution, and iteratively update the entries of $K$, respectively of $\Sigma$. Algorithms 1 and 2 are inspired by the corresponding algorithms for Gaussian graphical models; see, for example, [10, 39, 41]. Slawski and Hein [37] also provide a coordinate descent algorithm for estimating covariance matrices under ${\rm MTP}_{2}$. However, their method updates one column/row of $\Sigma$ at a time.

We first analyze Algorithm 1. Let $A=\{u,v\}$ and $B=V\setminus A$. Then note that the objective function can be written in terms of the $2\times 2$ Schur complement $K^{\prime}=K_{AA}-K_{AB}K_{BB}^{-1}K_{BA}$, since up to an additive constant

[TABLE]

Defining $L:=K_{AB}K_{BB}^{-1}K_{BA}$, then the optimization problem in step (2) of Algorithm 1 is equivalent to

[TABLE]

The unconstrained optimum to this problem is given by $K^{\prime}=S_{AA}^{-1}$ and is attained if and only if $(S_{AA}^{-1})_{12}+L_{12}\leq 0$, or equivalently, if and only if

[TABLE]

Otherwise the KKT conditions give that $K^{\prime}_{12}=-L_{12}$.

Maximizing over the remaining two entries of $K^{\prime}$ leads to a quadratic equation, which has one feasible solution

[TABLE]

Then the solution to the optimization problem in step (2) is given by $K_{AA}=K^{\prime}+L$.

Dual to this algorithm, one can define an algorithm that iteratively updates the off-diagonal entries of $\Sigma$ by maximizing the log-likelihood in direction $\Sigma_{uv}$ and keeping all other entries fixed. This procedure is shown in Algorithm 2. If $p>n$, $S$ is not positive definite; in this case we use as starting point the single linkage matrix $Z$ that is defined later in (13).

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\setminus A$ and $L=\Sigma_{AB}\Sigma_{BB}^{-1}\Sigma_{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

[TABLE]

We end by proving that Algorithms 1 and 2 indeed converge to the MLE. We here assume that $n\geq 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 $\hat{K}=\hat{\Sigma}^{-1}\in\mathcal{M}$.

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., [28]. 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 [26] 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 $K_{n}$ converging to $K$. Denote by $\tilde{K}$ and $\tilde{K}_{n}$ the corresponding one-step updates. We want to show that $\tilde{K}_{n}$ also converges to $\tilde{K}$. As above, let $A=\{u,v\}$, $B=V\setminus A$, $K^{\prime}=K_{AA}-K_{AB}K_{BB}^{-1}K_{BA}$ and $L=K_{AB}K_{BB}^{-1}K_{BA}$. Outside of the block $\tilde{K}_{AA}$ this convergence is trivial; so we focus only on the three entries in $\tilde{K}_{AA}$. The function $L_{12}\mapsto(K^{\prime}_{11},K_{22}^{\prime},K_{12}^{\prime})$ is continuous if and only if each coordinate is. It is clear that these functions are continuous if $L_{12}\neq\frac{S_{uv}}{S_{uu}S_{vv}-S_{uv}^{2}}$. It remains to show that if $L_{12}=\frac{S_{uv}}{S_{uu}S_{vv}-S_{uv}^{2}}$ the update in (11) gives $K^{\prime}=S_{AA}^{-1}$, which can be easily checked. ∎

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)

$U_{ii}\geq U_{ij}$ for all $i,j\in V$,

• (ii)

$U_{ij}\geq\min\{U_{ik},U_{jk}\}$ for all $i,j,k\in V$.

We say that a symmetric matrix is an inverse M-matrix if its inverse is an M-matrix. The connection between ultrametrics and 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$ is an inverse M-matrix.

The main reason why ultrametric matrices are relevant here is the following construction, which is similar to constructions used in in phylogenetics [34, Section 7.2] and single linkage clustering [16].

Let $R$ be a symmetric $p\times p$ positive semidefinite matrix such that $R_{ii}=1$ for all $i\in V$. Consider the weighted graph $G^{+}=G^{+}(R)$ over $V$ with an edge between $i$ and $j$ whenever $R_{ij}$ is positive and assign to each edge the corresponding positive weight $R_{ij}$. Note that $G^{+}$ in general does not have to be connected. Define a $p\times p$ matrix $Z$ by setting $Z_{ii}=1$ for all $i\in V$ and

[TABLE]

for all $i\neq j$, 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$.

Example 3.2**.**

Suppose that

[TABLE]

Then $G^{+}$ and $Z$ are given by

12340.60.50.40.2

[TABLE]

For example, to get $Z_{12}$ we consider two paths $1-3-2$ and $1-4-3-2$. The minimum of $R_{uv}$ over the first path is $0.4$ and over the second path $0.2$. This gives $Z_{12}=0.4$. ∎

Note that in the above example $Z\geq 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\times p$ positive semidefinite matrix satisfying $R_{ii}=1$ for all $i\in V$. Then the single-linkage matrix $Z$ based on $R$ is an ultrametric matrix with $Z_{ij}\geq R_{ij}$ for all $i\neq j$. If, in addition, $R_{ij}<1$ for all $i\neq 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, $R_{ij}\leq 1$ for all $i,j$ and from (13) it immediately follows that $Z_{ij}\leq 1$ and therefore $Z_{ii}\geq Z_{ij}$ for all $i,j$ as needed. Finally, to prove condition (ii) in the definition of ultrametric, let $i,j,k\in V$. Suppose first that $i,j,k$ lie in the same connected component of $G^{+}$. Let $P_{1}$, $P_{2}$ be the paths in $G^{+}$ such that $Z_{ik}=\min_{uv\in P_{1}}R_{uv}$ and $Z_{jk}=\min_{uv\in P_{2}}R_{uv}$. Let $P_{12}$ be the path between $i$ and $j$ obtained by concatenating $P_{1}$ and $P_{2}$. Then

[TABLE]

Now suppose that $i,j,k$ are not in the same connected component of $G^{+}$. In that case $0\in\{Z_{ij},Z_{ik},Z_{jk}\}$. Because zero is attained at least twice, again $Z_{ij}\geq\min\{Z_{ik},Z_{jk}\}$. Hence, $Z$ is an ultrametric matrix. The fact that $Z_{ij}\geq R_{ij}$ for all $i,j$ follows directly by noting that the edge $ij$ forms a path between $i$ and $j$.

Suppose now that $R_{ij}<1$ for all $i\neq j$. In that case also $Z_{ij}<1$ for all $i\neq 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\neq j$, then necessarily $Z_{ij}=Z_{ii}=Z_{jj}$, a contradiction. From Theorem 3.1 it then follows that $Z$ is an inverse M-matrix, which completes the proof. ∎

As a direct consequence we obtain the following result.

Proposition 3.4**.**

Let $S$ be a symmetric positive semidefinite matrix with strictly positive entries on the diagonal and such that $S_{ij}<\sqrt{S_{ii}S_{jj}}$ for all $i\neq j$. Then there exists an inverse M-matrix $Z$ such that $Z\geq S$ and $Z_{ii}=S_{ii}$ for all $i\in V$.

Proof.

Apply Proposition 3.3 to the normalized version $R$ of $S$, with entries $R_{ij}:=S_{ij}/\sqrt{S_{ii}S_{jj}}$. Because $R_{ij}<1$ for all $i\neq j$, the corresponding single-linkage matrix $Z^{\prime}$ is ultrametric with $Z^{\prime}\geq R$ and $Z^{\prime}$ is an inverse M-matrix. Define $Z$ by $Z_{ij}=\sqrt{S_{ii}S_{jj}}Z^{\prime}$. Then $Z\geq S$ and $Z_{ii}=S_{ii}$ for all $i\in V$. Moreover, $Z$ is an inverse M-matrix because $Z^{\prime}$ is. ∎

Proposition 3.4 is very important for our considerations. A basic application is an elegant alternative proof of the main result of [37], which says that the MLE under ${\rm MTP}_{2}$ exists with probability one as long as $n\geq 2$. This is in high contrast with the existence of the MLE in Gaussian graphical models without additional constraints; see [40].

Theorem 3.5** (Slawski and Hein [37]).**

Consider a Gaussian ${\rm MTP}_{2}$ model and let $S$ be the sample covariance matrix. If $S_{ij}<\sqrt{S_{ii}S_{jj}}$ for all $i\neq j$ then the MLE $\hat{\Sigma}$ (and $\hat{K}$) exists and it is unique. In particular, if the number $n$ of observations satisfies $n\geq 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\geq S$ and $Z_{ii}=S_{ii}$ for all $i$. It follows that $Z$ satisfies primal feasibility (6) and dual feasibility (7) and (8). By Proposition 2.2 the MLE exists and it is unique by convexity of the problem. ∎

Remark 3.6*.*

Combining this result with Corollary 2.4 we note that the cliques of $\hat{G}$ can at most be of size $n$. In this way the sparsity of $\hat{G}$ automatically adjusts to the sample size.

The matrix $Z$ can be computed efficiently222In our computations we use the single-linkage clustering method in R.. 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

[TABLE]

where $\overline{ij}$ denotes the unique path between $i$ and $j$ in $T$. For example $Z_{12}=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\times p$ positive semidefinite matrix satisfying $R_{ii}=1$ for all $i\in V$. Let ${\rm MWSF}(R)$ be the set of all minimal maximum weight spanning forests of $R$. Note that all edge weights of any such forest $F\in{\rm MWSF}(R)$ must be positive; hence we must have $F\subseteq G^{+}$. Also, if $R$ is an empirical correlation matrix, then ${\rm MWSF}(R)$ will be a singleton with probability one and in such cases we shall mostly speak of the MWSF.

Proposition 3.7**.**

The single-linkage matrix $Z$ as defined in (13) is block diagonal with blocks corresponding to the connected components of any $F\in{\rm MWSF}(R)$. Within each block all elements are strictly positive and given by

[TABLE]

where $\overline{ij}$ is the unique path between $i$ and $j$ in a maximal weight spanning tree of $R$. In particular, $Z_{ij}=R_{ij}$ for all edges of ${\rm MWSF}(R)$.

Proof.

First suppose that $i,j\in V$ lie in two different components of $F\in{\rm MWSF}(R)$. This means that there is no path between $i$ and $j$ in $G^{+}$ and so, by definition, $Z_{ij}=0$. Because $Z_{ij}>0$ if $i,j$ lie in the same component of $F$, $Z$ is block diagonal with blocks corresponding to connected components of ${\rm MWSF}(R)$.

The rest of the proof is an adaptation of a proof of a related result [34, Proposition 7.2.10]. Suppose that $i,j\in V$ lie in the same connected component of $F$ and denote the tree in $F$ corresponding to this component by $T$. By definition $Z_{ij}\geq\min_{uv\in\overline{ij}}R_{uv}$. Suppose that $Z_{ij}>\min_{uv\in\overline{ij}}R_{uv}$. 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 $Z_{ij}>R_{kl}$, there exists a path $P$ in $G^{+}$ between $i$ and $j$ such that $R_{uv}>R_{kl}$ 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^{\prime}l^{\prime}$ (distinct from $kl$) in $P$ whose end vertices lie in different blocks of this partition. Let $T^{\prime}$ be the spanning tree obtained from $T$ by deleting $kl$ and adding $k^{\prime}l^{\prime}$. Since $R_{k^{\prime}l^{\prime}}>R_{kl}$, the total weight of $T^{\prime}$ is greater than $T$, which is a contradiction. We conclude that $Z_{ij}=\min_{uv\in\overline{ij}}R_{uv}$ for all $i,j$ in the same connected component of $G^{+}$. ∎

To conclude this section, we note that the starting point $\Sigma^{0}$ of Algorithm 2 is arbitrary as long as $\Sigma^{0}\succ 0$ and $\Sigma^{0}\geq S$. The single-linkage matrix $Z$ 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\notin\mathbb{S}_{\succ 0}$ and hence not feasible.

4 The maximum likelihood graph

Fitting a Gaussian model with ${\rm MTP}_{2}$ constraints tends to induce sparsity in the maximum likelihood estimate $\hat{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 $R_{ii}=1$ for all $i$ and $R_{ij}<1$ for all $i\neq j$. Consider again the weighted graph $G^{+}=G^{+}(R)$. 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 $\hat{\Sigma}$ under ${\rm MTP}_{2}$ is a block diagonal matrix with strictly positive entries in each block. The blocks correspond precisely to trees in ${\rm MWSF}(R)$.

Proof.

Firstly, since $\hat{\Sigma}$ is an inverse M-matrix, it is block diagonal with strictly positive entries in each block; see, e.g., Theorem 4.8 in [22]. We will show that each block of $\hat{\Sigma}$ corresponds precisely to a tree in ${\rm MWSF}(R)$.

Denote the vertex sets for a forest $F\in{\rm MWSF}(R)$ as $T_{1},\ldots,T_{k}$ and the blocks of $\hat{\Sigma}$ as $B_{1},\ldots,B_{l}$. Firstly, for any $T_{i}$ there must be a $j$ so that $T_{i}\subseteq B_{j}$; this is true since all entries in $R$ along the edges of $T_{i}$ are positive and thus $\hat{\Sigma}\geq R>0$. Thus the block partitioning corresponding to the trees is necessarily finer than that of $\hat{\Sigma}$.

On the other hand, suppose that two different trees $T_{i}$ and $T_{j}$ in $F$ are in the same block of $\hat{\Sigma}$ so that $\hat{\Sigma}_{uv}>0$ for all $u\in T_{i}$ and $v\in T_{j}$. Then, as we must have $R_{uv}\leq 0$, also necessarily $\hat{\Sigma}_{uv}-R_{uv}>0$. Complementary slackness (9) now implies that $\hat{K}_{uv}=0$ for all $u\in T_{i}$ and $v\in T_{j}$, and hence $\hat{K}$ is block-diagonal with blocks corresponding to the trees in $F$. Since $\hat{\Sigma}=\hat{K}^{-1}$, we also get $\hat{\Sigma}_{uv}=0$ which contradicts that $u$ and $v$ are in the same block of $\hat{\Sigma}$. ∎

This result shows that, without loss of generality, we can always assume that $G^{+}$ is connected and then ${\rm MWSF}(R)={\rm 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 ${\rm MWSF}(R)$ are just trees.

4.1 An upper bound on the ML graph

In the following, we provide a simple procedure for identifying an upper bound for $\hat{G}$. This procedure relies on the estimation of the standard Gaussian graphical model over the tree ${\rm MWSF}(R)$. The MLE under this assumption, denoted by $\tilde{\Sigma}$, can be computed efficiently and it satisfies

[TABLE]

where $\overline{ij}$ denotes the unique path between $i$ and $j$ in ${\rm MWSF}(R)$; see, for example, [42, Section 8.2].

To provide an upper bound on $\hat{G}$, we will make use of a connection to so-called path product matrices: A non-negative matrix $R$ is a path product matrix if for any $i,j\in V$, $k\in\mathbb{N}$, and $1\leq i_{1},\ldots,i_{k}\leq p$

[TABLE]

If in addition the inequality is strict for $i=j$, we say that $R$ is a strict path product matrix. We note the following:

Theorem 4.2** (Theorem 3.1, [21]).**

Every inverse M-matrix is a strict path product matrix.

We are now able to provide an upper bound for the ML graph $\hat{G}$.

Proposition 4.3**.**

*The pair $ij$ forms an edge in the ML graph only if

[TABLE]

for any path in $P$ in $G^{+}$. In particular, $R_{ij}\leq 0$ implies that $ij$ is not an edge of the ML graph.

Proof.

Because $\hat{\Sigma}$ is an inverse M-matrix it is necessarily a path product matrix by Theorem 4.2. In particular, for all $i,j$

[TABLE]

for any path $P$ in $G^{+}$. By Proposition 2.2, we also have $\hat{\Sigma}_{uv}\geq R_{uv}$. Thus, if $ij\in\hat{G}$

[TABLE]

as desired. ∎

Motivated by this result we define the excess correlation graph (EC graph) ${\rm EC}(R)$ of $R$ by the condition

[TABLE]

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,

[TABLE]

The inclusion is typically strict. For example, if $R$ is an inverse M-matrix, then ${\rm EC}(R)$ is the complete graph, whereas $G(\hat{K})$ can be arbitrary; this follows from [13, Proposition 6.3].

4.2 Some exact results on the ML graph

Next, we analyze generalization of trees known as block graphs, where edges are replaced by cliques, and give a condition under which the maximum likelihood estimator admits a simple closed-form solution. More formally, $G$ is a block graph if it is a chordal graph with only singleton separators. It is natural to study block graphs, since viewing the MLE $\hat{\Sigma}$ as a completion of $S$, block graphs play the same role for inverse M-matrices as chordal graphs play for Gaussian graphical models, see for example [20] and Corollary 7.3 of [13].

We first define a matrix $W=W(R)$ by

[TABLE]

where, like in (13), the maximum is taken over all paths in $G^{+}$ between $i$ and $j$ and is set to zero if no such path exists. Transforming $D_{ij}=-\log R_{ij}$ gives a distance based interpretation, in which $W_{ij}$ is related to the shortest distance between $i$ and $j$ in $G^{+}$ with edge lengths given by $D_{uv}$. We also have the following simple lemma.

Lemma 4.4**.**

The matrix $W$ is a path product matrix. Further, $R$ is a path product matrix if and only if $W(R)=R$.

Proof.

This is immediate from the definition of $W$.∎

It is easy to show that $Z\geq W\geq R$ and that $W$ is always equal to the MLE $\hat{\Sigma}$ in the case when $p\leq 3$. For general $p$ we do not know conditions on $R$ that assure that $W$ is an inverse M-matrix, or the MLE. Indeed, Example 3.4 in [21] gives a strict path product correlation matrix $R$, and thus $W=R$, which is not an inverse M-matrix, and thus $W\neq\hat{\Sigma}$. We note that $W=\hat{\Sigma}$ for the carcass data discussed in Example 1.1 and, as we shall see in the following, it reflects that in this example, the ML graph is a block graph.

Let $G_{R}(W)$ be the graph having edges $ij$ exactly when $R_{ij}=W_{ij}$ and no edges otherwise. We then obtain the following result.

Proposition 4.5**.**

If $G_{R}(W)$ is a block graph and blocks of $W$ corresponding to cliques are inverse M-matrices, then $\hat{\Sigma}=W$ and $\hat{G}\subseteq G_{R}(W)$.

Proof.

Note first that if $\hat{\Sigma}=W$, the KKT conditions (9) imply that $\hat{G}\subseteq G_{R}(W)$. Let $\tilde{\Sigma}$ denote the maximum likelihood estimate of $\Sigma$ under the Gaussian graphical model with graph $G_{R}(W)$. Then, since $G_{R}(W)$ is a block graph, it follows from [26, equation (5.46) on page 145] that $\tilde{\Sigma}$ is an inverse M-matrix which coincides with $W$ and $R$ on all edges of $G_{R}(W)$. So from to show that $\tilde{\Sigma}=\hat{\Sigma}=W$ we just need to argue that $\tilde{\Sigma}=W$.

We proceed by induction on the number $m=|\mathcal{C}|$ of cliques of $G_{R}(W)$. If there is only one clique in $G_{R}(W)$, we have $\tilde{\Sigma}=R$ and $R$ is an inverse M-matrix and hence $\hat{\Sigma}=R=W$. Assume now that the statement holds for $|\mathcal{C}|\leq m$ and assume $G_{R}(W)$ has $m+1$ cliques. Since $G_{R}(W)$ is a block graph, there is a decomposition $(A,B,S)$ of $G_{R}(W)$ into block graphs with at most $m$ cliques and with the separator $S=\{s\}$ being a singleton. But for a decomposition of $G_{R}(W)$ as above we have from [26, equation (5.31) in Proposition 5.6] and the inductive assumption that

[TABLE]

Now let $P^{*}$ be the path in $G^{+}$ such that $W_{ij}=\prod_{uv\in P^{*}}R_{uv}$ for any two vertices $i,j$. We claim that all edges in $P^{*}$ must be edges of $G_{R}(W)$. Otherwise, suppose $P^{*}$ contains an edge $uv$ which is not an edge in $G_{R}(W)$; then $W_{uv}>R_{uv}$ and so if we replace the edge $uv$ with the path realizing $W_{uv}$ the product would be strictly increased, which contradicts the optimality of $P^{*}$. Since $S$ is a singleton separator, this also implies that $P^{*}$ passes through $S$ whenever it involves vertices from both $A$ and $B$. Suppose that $i,j\in A\cup S$. Then optimality of $P^{*}$ implies that $P^{*}$ is contained in $A\cup S$ and so $\tilde{\Sigma}_{A\cup S}=W(R_{A\cup S})=W_{A\cup S}$ and by the same argument $\tilde{\Sigma}_{B\cup S}=W_{B\cup S}$. Moreover, if $i\in A$ and $j\in B$ then $W_{ij}=W_{is}W_{sj}$. Now the inductive assumption in combination with the expression [26, page 140] yields that

[TABLE]

and thus $\tilde{\Sigma}=\hat{\Sigma}=W$ as required. ∎

Remark 4.6*.*

We note that with probability one, the slackness constraints in (9) are not simultaneously active, and hence in Proposition 4.5 we have almost sure equality between $G_{R}(W)$ and $\hat{G}$. Thus we can identify $\hat{G}$ without first calculating $\hat{K}$.

We further have the following corollary.

Corollary 4.7**.**

Under the same conditions as in Proposition 4.5 we have ${\rm MWSF}(R)\subseteq\hat{G}\subseteq G_{R}(W)$.

Proof.

Consider an edge $ij$ between vertices in different cliques of $G_{R}(W)$ and assume $S_{1}=\{s_{1}\},S_{2}=\{s_{2}\}$ are $(i,j)$-separators with $i\sim s_{1}$ and $j\sim s_{2}$. Then, since $\hat{G}\subseteq G_{R}(W)$ we have i\mbox{\,\perp!!!\perp\,}j\,|\,s_{1} and i\mbox{\,\perp!!!\perp\,}j\,|\,s_{2} according to $\hat{\Sigma}$ and therefore

[TABLE]

so the edge $ij$ can never be part of a ${\rm MWSF}$ because removing the edge would render either $s_{1}$ disconnected from $i$ or $s_{2}$ disconnected from $j$ and then the weight of a ${\rm MWSF}$ would increase when replacing $ij$ with $is_{1}$ or $js_{2}$, respectively. This completes the proof. ∎

It is not correct in general that ${\rm MWSF}(R)\subseteq\hat{G}$ as demonstrated in the following example; although this has been the case in all non-constructed examples we have considered including the relatively large Example 5.8 below.

Example 4.8**.**

Below we display an inverse M-matrix $K$

[TABLE]

and the corresponding correlation matrix

[TABLE]

Here the ${\rm MWSF}(R)$ is the star graph with $5$ as its center, but the edge $2\sim 5$ is not in $G(K)$. Note that all the edges in $G^{+}$ adjacent to $2$ have almost the same weight. We note that we have also calculated $K^{-1}$ using rational arithmetic to ensure the phenomenon cannot be explained by rounding error.

5 Gaussian signed ${\rm MTP}_{2}$ distributions

In this section we discuss how our results can be generalized to so-called signed ${\rm MTP}_{2}$ Gaussian distributions, where the distribution is ${\rm MTP}_{2}$ up to sign swapping. Such distributions were discussed by Karlin and Rinott [24]. More precisely, a random variable $X$ has a signed ${\rm MTP}_{2}$ distribution if there exists a diagonal matrix $D$ with $D_{ii}=\pm 1$ (called sign matrix) such that $DX$ is ${\rm MTP}_{2}$. The following characterization of signed ${\rm MTP}_{2}$ Gaussian distributions is a direct consequence of [24, Theorem 3.1 and Remark 1.3].

Proposition 5.1**.**

A Gaussian random variable $X$ has a signed ${\rm MTP}_{2}$ distribution if and only if $|X|$ is ${\rm MTP}_{2}$.

Gaussian graphical models with signed ${\rm MTP}_{2}$ distributions are called non-frustrated in the machine learning community. The following result is implicitly stated in [29].

Theorem 5.2**.**

A Gaussian random variable $X$ with concentration matrix $K$ has a signed ${\rm MTP}_{2}$ distribution if and only if it holds for every cycle $(i_{1},\ldots,i_{k},i_{1})$ in the graph $G(K)$ that

[TABLE]

Proof.

The ’only if’ direction is easy to check. Note that (15) can be rephrased by saying that each cycle in the graph with edge weights given by the off-diagonal entries of $-K$ has an even number of negative edges. The ’if’ direction can now be recovered from the proof of [29, Corollary 3]. ∎

Signed ${\rm MTP}_{2}$ 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 ${\rm MTP}_{2}$ distributions. The ${\rm MTP}_{2}$ distributions among those are precisely those without negative entries in the covariance matrix $\Sigma$.

Proof.

Let $T$ be a tree and $K=\Sigma^{-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 ${\rm MTP}_{2}$. The second part of the statement follows from [13, Corollary 7.3]. ∎

Because signed ${\rm MTP}_{2}$ 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 [24].

Theorem 5.4**.**

Every distribution on a Gaussian tree model with hidden variables is signed ${\rm MTP}_{2}$.

Gaussian tree models with hidden variables have many applications, in particular related to modeling evolutionary processes; see, e.g., [7, 36]. As an important submodel they contain the Brownian motion tree model [14]. 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 ${\rm MTP}_{2}$ distributions in this model correspond to the distributions in a Spearman model [27, 38], where the hidden factor is interpreted as intelligence.

Let $R$ be a sample correlation matrix. Maximizing the likelihood over all signed ${\rm MTP}_{2}$ Gaussian distributions requires determining the sign matrix $D$, with $D_{ii}=\pm 1$, that maximizes the likelihood for all $2^{p}$ possible matrices $DRD$. A natural heuristic is to choose $D^{*}$ such that $D_{ii}^{*}D_{jj}^{*}R_{ij}\geq 0$ for all edges $ij$ of ${\rm MWSF}(|R|)$, where $|R|$ denotes the matrix whose entries are the absolute values of the entries of $R$. We provide conditions under which this procedure indeed leads to the MLE under signed ${\rm MTP}_{2}$, and we also provide examples showing that this is not true in general. Quite interestingly, balanced graphs again play an important role in this part of the theory.

First we describe how to obtain a sign swapping matrix $D^{*}$ such that $D_{ii}^{*}D_{jj}^{*}R_{ij}\geq 0$ for all edges $ij$ of ${\rm MWSF}(|R|)$. Root ${\rm MWSF}(|R|)$ at node $1$, that is, regard ${\rm MWSF}(|R|)$ as a directed tree with all edges directed away from $1$. Set $D_{11}^{*}=1$. Then proceed recursively. For any edge $i\to j$ suppose that $D_{ii}^{*}$ is known and set $D_{jj}^{*}:={\rm sgn}(D^{*}_{ii}R_{ij})$. Note that by construction

[TABLE]

where $1\to i_{1}\to i_{2}\to\cdots\to i_{k}\to i$ is the unique path from $1$ to $i$ in ${\rm MWSF}(|R|)$. We set $D^{*}_{ii}=0$ if no such path exists. It is easy to check that the resulting $D^{*}$ satisfies $D_{ii}^{*}D_{jj}^{*}R_{ij}\geq 0$ for all edges $ij$ of ${\rm MWSF}(|R|)$.

Proposition 5.5**.**

Suppose that $R$ is a sample correlation matrix whose graph is balanced, that is, such that for every cycle $(i_{1},i_{2},\ldots,i_{k},i_{1})$ in the graph $G(R)$

[TABLE]

Then the MLE based on $R$ over signed ${\rm MTP}_{2}$ Gaussian distributions is equal to the MLE based on the sample correlation matrix $D^{*}RD^{*}$ over ${\rm MTP}_{2}$ distributions.

Proof.

We first show that $D^{*}RD^{*}$ has only positive entries. Let $i,j$ be any two nodes and let $1\to i_{1}\to\cdots\to i_{k}\to i$ and $1\to j_{1}\to\cdots\to j_{l}\to j$ be the paths in ${\rm MWSF}(|R|)$ from $1$ to $i$ and $j$, respectively. By (16) we obtain

[TABLE]

which is positive by (17). This shows that without loss of generality we can assume that all entries of $R$ are nonnegative and hence that $D^{*}$ is the identity matrix $\mathbb{I}_{p}$. We now show that the likelihood over ${\rm MTP}_{2}$ distributions given the sample correlation matrix $DRD$ is maximized by $D=\mathbb{I}_{p}$. This is because $(D_{ii}D_{jj}-1)\leq 0$ and $R_{ij}K_{ij}\leq 0$, and hence

[TABLE]

which completes the proof. ∎

Note that any spanning tree $T$ of $G^{+}(|R|)$ would suffice to identify the sign switches as above.

Proposition 5.5 provides a sufficient condition for $D^{*}$ to be the optimal sign-switching matrix; i.e., it provides a sufficient condition such that for every $K\in\mathbb{S}_{\succ 0}$ and every sign matrix $D$ it holds that

[TABLE]

As a consequence of Proposition 5.5 we obtain the following result for the case when the sample size is 2.

Corollary 5.6**.**

If the sample correlation matrix $R$ is based on $n=2$ observations, then the MLE over signed ${\rm MTP}_{2}$ Gaussian distributions given $R$ is equal to the MLE over ${\rm MTP}_{2}$ Gaussian distributions given the modified sample correlation matrix $D^{*}RD^{*}$.

Note that the case $n=2$ is special and Proposition 5.5 does not extend to arbitrary sample correlation matrices. In the following, we give a simple counterexample.

Example 5.7**.**

Suppose that the sample correlation matrix is

[TABLE]

Then ${\rm MWSF}(|R|)$ is given by the star graph with edges $1-2$, $1-3$, $1-4$. Since $R$ is positive on these entries, $D^{*}=\mathbb{I}_{p}$. 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

[TABLE]

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

[TABLE]

Note that $R_{13}$ is only slightly larger than $-R_{23}$ and $-R_{24}$. Hence, in essence we are increasing the number of constraints by one, which explains the decrease of the log-likelihood value. ∎

We conclude 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 [30]333We downloaded the data from http://web.stanford.edu/class/psych253/tutorials/FactorAnalysis.html.: 240 individuals were asked to rate themselves on the scale 1-9 with respect to 32 different personality traits. The resulting correlation matrix is shown in Figure 2.

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 set follows a signed ${\rm MTP}_{2}$ distribution and analyze it under this constraint.

The correlation matrix resulting from the sign switching procedure described in (16) is shown on the left in Figure 3, while the correlation matrix resulting from switching the signs of the 16 (negative) traits that constitute the first block of variables in Figure 2 is shown on the right in Figure 3.

These plots suggest that the matrix on the right is closer to being ${\rm MTP}_{2}$. In fact, its log-likelihood (i.e., the value of $\,\frac{n}{2}(\log\det K-{\rm tr}(SK))$) is -2046.146, as compared to the log-likelihood value of -2071.717 resulting from the sign switching procedure described in (16). For comparison, the value of the unconstrained log-likelihood is -1725.075 and the value of the log-likelihood under ${\rm MTP}_{2}$ without sign switching is -2356.639. The unconstrained log-likelihood gives a lower bound of 642.142 on the likelihood ratio statistic to test signed ${\rm MTP}_{2}$ constraints, while the likelihood ratio statistic to test ${\rm MTP}_{2}$ constraints against the saturated model is equal to 1263.128.

The graphical models based on no sign switching and switching the signs of the 16 negative traits are shown in Figure 4.

The vertex labels are as shown in Table 1.

The red edges correspond to the maximum weight spanning trees. Red and blue edges together form the edge set of the ML graph so in both of these cases we have ${\rm MWSF}(R)\subseteq\hat{G}$. Finally, the grey edges are the remaining edges in the EC graph. As expected, the graph on the right looks denser. The interpretation of the spanning tree in both cases is very different. Edges in the first one connect similar personalities such as 6-24 (agreeable and cooperative), 12-22 (outgoing and sociable), 11-23 (disorganized and lazy). On the other hand, the second tree looks similar but it links also some almost perfect opposite personalities such as 12-14 (outgoing and shy), 22-30 (sociable and withdrawn), 11-26 (disorganized and organized), 7-10 (tense and relaxed). Note that none of these four edges are part of the ML graph on the left in Figure 4.

6 Discussion

In this article we have investigated maximum likelihood estimation for Gaussian distributions under the restriction of multivariate total positivity, used a connection to ultrametrics to show that it has a unique solution when the number of observations is at least two, shown that under certain circumstances the MLE can be obtained explicitely, and given convergent algorithms for calculating the MLE. For signed ${\rm MTP}_{2}$ distributions we have also given conditions under which a heuristic procedure for applying sign changes is correct and can be used to obtain the MLE.

It remains an issue to consider the asymptotic properties of the estimators we have given, and to derive reliable methods for identifying whether a given sample is consistent with the ${\rm MTP}_{2}$ assumption. On the former issue, standard arguments for convex exponential families ensure that if the true value $K_{0}$ is an M-matrix, $\hat{K}$ is a consistent estimator of $K_{0}$; and this is true whether or not the ${\rm MTP}_{2}$ assumption is envoked.

Another question is whether the ML graph $\hat{G}$ will be consistent for the true dependence graph. It is clear that without some form of penalty or thresholding, it cannot be the case. For example, if $p=2$ and the true $\Sigma$ is a diagonal matrix, the distribution of the empirical correlation $R_{12}$ will be symmetric around [math]. Hence, with probability $1/2$ the ML graph contains an edge between $1$ and $2$ and with probability $1/2$ it does not contain such an edge. This phenomenon persists for any number of observations $n$. Thus, to achieve consistent estimation of the dependence graph of $\Sigma$, some form of penalty for complexity or thresholding must be applied, the latter being suggested by [37], who also suggest a refitting after thresholding to ensure positive definiteness of the thresholded matrix. However, positive definiteness is automatically ensured, as shown below.

Proposition 6.1**.**

Let $K$ be an M-matrix over $V$ and $G=(V,E)$ an undirected graph. Define $K^{G}$ by

[TABLE]

Then $K^{G}$ is an M-matrix.

Proof.

We may without loss of generality assume that $K$ is scaled such that all diagonal elements are equal to 1; also it is clearly sufficient to consider the case when only a single off-diagonal entry $K_{uv}$ is replaced by zero. We have to show that the resulting matrix $K^{G}$ is positive definite.

Now, let $A=\{u,v\}$ and $B=V\setminus A$ and consider the Schur complements

[TABLE]

Since $K^{G}_{BB}=K_{BB}$, $K^{G}$ is positive definite if and only if $K^{G}/K_{BB}$ is. Because $K$ is an M-matrix, all entries in $K_{AB}(K_{BB})^{-1}K_{BA}$ are non-negative. Hence, we can write the Schur complements as

[TABLE]

where $c,d\in(0,1)$ and $a,b\geq 0$. Since $K$ is positive definite we have

[TABLE]

and hence

[TABLE]

implying that $K^{G}/K_{BB}$ is positive definite. This completes the proof. ∎

The consistency of the estimator $\hat{K}$ ensures that the ML graph will eventually contain the true dependence graph when $n$ becomes large and with an appropriate thresholding or penalization, this ensures that the true graph can be recovered, as also argued in [37].

The issue of the asymptotic distribution of the likelihood ratio test for ${\rm MTP}_{2}$ is an instance of testing a convex hypothesis within an exponential family of distributions. In our particular case, the convex hypothesis is a polyhedral cone with facets determined by the dependence graph $G(K)$. In such cases, the likelihood ratio test for the convex hypothesis typically has an asymptotic distribution which is a mixture of $\chi^{2}$-distributions with degrees of freedom determined by the co-dimension of these facets; see for example the analysis of the case of multivariate positivity in models for binary data by [3], using results of [35].

While these issues are both interesting and important, we consider them to be outside the scope of the present paper as they may be most efficiently dealt with in the more general context of exponential families, containing both the Gaussian and binary cases as special instances. We plan to return to these and other problems in the future.

Acknowledgements

Caroline Uhler was partially supported by DARPA (W911NF-16-1-0551), NSF (DMS-1651995), ONR (N00014-17-1-2147), and a Sloan Fellowship. We also thank two anonymous referees for their helpful comments.