Multivariate Generalized Gaussian Distribution: Convexity and Graphical Models

Covariance estimation is a fundamental problem in multivariate statistics. Many techniques for hypothesis testing, inference, denoising and prediction rely on accurate estimation of the true covariance. The problem is challenging when the available data is high dimensional and non-Gaussian. Such settings are typical in many applications including speech, radar, wireless communication, finance and more. These led to a growing interest in both robust and structured covariance estimation. Specifically, in this paper, we consider maximum likelihood estimation (MLE) in the multivariate generalized Gaussian distribution, with and without sparsity constraints on the inverse covariance.

The first part of this paper considers the geodesic convexity in MLE in elliptically symmetric (ES) distributions. Methods for robust covariance estimation date back to the early works of [18, 31]. A popular approach is Tyler’s scatter estimate. It involves a non-convex optimization yet can be solved via simple fixed point iteration. It has been rigorously analyzed and successfully applied to different problems [31, 26]. Recently, it was shown that the result is in fact the solution to a geodesically convex minimization [6, 10, 32, 33, 36]. Geodesic convexity ensures that any local minima is also globally optimal and leads to a much simpler analysis [27]. It also allows for numerous extensions, e.g., regularized solutions. In a competing line of works, a different class of robust covariance estimation techniques was proposed based on the MGGD [23, 11]. A well known example of MGGD is the multivariate Laplace distribution [14]. Fixed point iterations for MGGD estimation and their analyses has recently been considered in [24, 9, 25]. The first contribution in this paper is a new analysis which shows that the negative log-likelihood in MGGD is also geodesically convex. This result requires weaker conditions than previous analyses, provides more intuition and paves the road to numerous generalizations.

The second part of this paper addresses structured covariance estimation. Structure exploitation is a main ingredient in modern statistics that allows accurate high dimensional estimation via a small number of samples. A promising approach is based on sparse inverse covariance models. In the multivariate Gaussian case these are known as graphical models and characterize conditional independence [13, 20, 8, 4]. These models have been successfully applied to speech recognition, sensor networks, computer networks and other fields in signal processing [10, 35]. Our goal is to combine such models with non-Gaussian distributions, e.g., MGGD. Recently, a similar problem was addressed using an expectation maximization technique [15]. Another line of work focused on combining Tyler’s scatter estimate with a banded inverse covariance prior [2, 1]. Graphical models for transelliptical distributions were discussed in [21]. In this paper, we combine the MGGD framework with prior sparsity constraints on the inverse covariance. We show that the optimization can be formulated into a convex form as long as the MGGD scale parameter is larger than half and the sparsity satisfies a chordal structure. Chordal models, also known as decomposable or triangulated models, include banded structures, multiscale settings and other practical scenarios [34, 16, 12]. Such structures are associated with a perfect ordering of the variables. A typical example is banded models associated with time-varying autoregressive processes [3].

Our results are also applicable to the case of unknown sparsity pattern, i.e., structure learning via sparsity inducing penalties, but require prior knowledge of the perfect order. Recent works on structure learning in directed acyclic graphs provide data driven techniques for learning this order [30, 28].

The outline of the paper is as follows. We begin in Section II with a few mathematical preliminaries that will be useful in our work. Then, we continue to our two main contributions. In Section III we provide a new geodesic analysis of MGGD estimation, and in Section IV we introduce a convex optimization framework for chordal structured MGGD estimation. Simulation results are described in Section V, and concluding remarks are offered in Section VI.

We use the following notations. We denote the set of real, symmetric and positive definite matrices by $\mathrm{S_{++}}(p)\subset\mathbb{R}^{p\times p}$. We denote the span operator by $\mathrm{sp}\{\cdot\}$.

We begin with a brief review of two mathematical concepts which will be instrumental in the next sections.

In this section, we consider unconstrained MLE in MGGD. More precisely, we address a more general family of elliptically symmetric (ES) distributions (see [11] for a review and [24] for a recent generalization to complex case). These problems involve non-convex minimizations, yet it has been shown that their global minima can be efficiently found using simple fixed point iterations. We will show that the negative log-likelihoods are in fact geodesically convex, and that this may be the underlying principle behind their success.

An random variable $\bm{z}\in\mathbb{R}^{p}$ has a ES distribution in real space if its probability density function (p.d.f.) is

where $g:[0,\infty)\rightarrow(0,\infty)$ such that $\int_{0}^{\infty}t^{p-1}g(t^{2}){\,\mathrm{d}}t<\infty$, $C_{p,g}$ is a normalization constant such that the integral of the distribution is $1$. In (5), $\bm{\Sigma}\in\mathrm{S_{++}}(p)$ is called a scatter matrix, and $\mu\in\mathbb{R}^{p}$ is the center of the distribution.

MGGD [23] is a widely used special case of ES when

where $\beta$ is the shape parameter. In particular, for $\beta=0.5$ it gives multivariate analog of Laplace distribution, and the multivariate Gaussian distribution is obtained for $\beta=1$.

We consider the estimation of $\bm{\Sigma}$ given $n$ independent and identically distributed (i.i.d) realizations of a zero mean ES random vector denoted by $\bm{z}_{1},\cdots,\bm{z}_{n}$. The MLE is the parameter that minimizes the negative-log-likelihood

where $\rho(x)=-\log g(x)$. The following Theorem 1 characterizes the existence and uniqueness properties of this MLE. Theorem 1(a) characterizes the uniqueness and is the main contribution in this section; Theorem 1(b) characterize the existence and is borrowed from [31, Theorem 2.1].

Before proving this theorem, a few comments are in order. We remark that the condition $\mathrm{sp}\{\bm{z}_{1},\bm{z}_{2},\cdots,\bm{z}_{n}\}=\mathbb{R}^{p}$ is also implicated by (7), since otherwise $\mathrm{V}=\mathrm{sp}\{\bm{z}_{1},\bm{z}_{2},\cdots,\bm{z}_{n}\}$ violated the assumption. This condition is necessary since otherwise $L_{0}(\bm{P}_{\mathrm{V}})\rightarrow-\infty$ as $\bm{\Sigma}\rightarrow\bm{P}_{\mathrm{V}}$, where $\bm{P}_{\mathrm{V}}$ is the projector to subspace $\mathrm{V}$) and therefore the minimizer of $L_{0}(\bm{\Sigma})$ does not exist.

Theorem 1 relaxes the conditions for the uniqueness/existence of the minimizer of $L_{0}(\bm{\Sigma})$ presented in [19, Theorem 2.2]:
M1 $\rho^{\prime}(x)$ is non-negative, continuous and nonincreasing.
M2 $x\rho^{\prime}(x)$ is strictly increasing.
M3 condition (7).
Our assumption that $\rho$ is increasing corresponds to $\rho^{\prime}(x)>0$, which follows from M1, where $\rho^{\prime}$ is nonnegative and M2 (which exclude the possibility that $\rho^{\prime}(x)=0$ for $x>0$); and our assumption that $\rho(e^{x})$ is strictly convex corresponds to M2, where $x\rho^{\prime}(x)$ is strictly increasing. In comparison, our conditions does not require $\rho$ to be differentiable, and we do not assume their assumptions in M1 that $\rho^{\prime}(x)$ is continuous or nonincreasing.

A special case of Theorem 1 is Tyler’s M-estimator in which $\rho(x)=\log(x)/2p$. Following the first part of Theorem 1(a), $L_{0}(\bm{\Sigma})$ is geodesically convex, which has been previously identified in [6, 32, 33, 36]. The geodesic convexity does not contradict the non-uniqueness of Tyler’s M-estimator since this convexity is not strict. Another special case is the class of MGGD estimators. The following corollary follows from Theorem 1 and the fact that $a_{1}$ in Theorem 1(b) is $\infty$.

Existence and uniqueness of the MLE in MGGD has been previously addressed in [24, Section V.A]. However, this contribution does not identify geodesic convexity and applies only to $0<\beta\leq 1$. Our theorem applies to MGGD with all $\beta>0$ and provides additional insight based on geodesic convexity. Furthermore, it applies to other ES distributions, including cases where $\rho^{\prime}(x)$ is not continuous, e.g., $\rho(x)=x$ when $x\leq 1$ and $\rho(x)=3x-2$ when $x>1$.

Here we remark that although ML estimator for elliptically symmetric (ES) distribution is the motivation of the argument, Theorem 1 (and Theorem 2 in next section) hold even if $\rho(x)$ does not relate to an elliptically symmetric (ES) distribution. For example, Tyler’s M-estimator can be written in the form of (5) with $\rho(x)=\log(x)/2p$; however this $\rho(x)$ corresponds to the central angular Gaussian which is not a member in ES class [24, (36)]. Another example is in [24, page 5610, Example 1], where $\rho^{\prime}$ is set to be a huber function.

We remark that the Theorem 1 also applied to the complex elliptically symmetric (CES) distributions defined by

where the MLE estimator minimizes

Then Theorem 1 holds with $a_{1}=\sup\{a|x^{a}\exp(-\rho(x))\rightarrow 0\,\,\,\,\text{as $x\rightarrow\infty$}\}$. This is also a generalization of the uniqueness/existence of the minimizer of $L_{0}(\bm{\Sigma})$ presented in [24, Section V.A].

In practice, various numerical techniques can be used to find a local minima of the MGGD negative log-likelihood. Theorem 1 and geodesic convexity then ensure that this local minima will also be global minima. A promising approach is the classical iterative reweighed scheme due to [19, 5]:

where $u(x)=\rho^{\prime}(x)$. It has been shown in [5, Proposition 1] that when $\mathrm{Sp}(\bm{z}_{1},\bm{z}_{2},\cdots,\bm{z}_{n})=\mathbb{R}^{p}$, $\rho^{\prime\prime}(x)$ is continuous and nonnegative, then $L_{0}(\bm{\Sigma}_{m})$ decreases monotonically ($\mathrm{Sp}(\bm{z}_{1},\bm{z}_{2},\cdots,\bm{z}_{n})=\mathbb{R}^{p}$ is assumed in order to make sure that $L_{0}(\bm{\Sigma}_{m+1})$ is well-defined). [5, Proposition 3] shows that any limiting point of the sequence $\bm{\Sigma}_{m}$ is a stationary point. When the assumptions in Theorem 1 hold, this point is the unique minimizer of $L_{0}(\bm{\Sigma})$.

In this section, we consider structured MGGD estimation. In particular, we consider MLE of the MGGD scatter matrix subject to sparsity constraints. Thus, we are interested in the solution to

where the objective is the negative MGGD log-likelihood as described in Section III, and $E$ is the edge set of a known graph $G(V,E)$ associated with the sparsity of $\bm{\Sigma}^{-1}$. Unfortunately, the above minimization is not convex in $\bm{\Sigma}$ or $\bm{\Sigma}^{-1}$. Nor is it geodesically convex in it. Indeed, the sparsity constraints are not preserved by the positive definite geodesic in (3). Thus, it is not clear whether its global minima can be found in an efficient manner.

In what follows, we propose a simple trick to “convexify” the optimization in many interesting cases of MGGD. In particular, we assume the $G$ is chordal and represent $\bm{\Sigma}^{-1}$ with its Cholesky factor $\bm{C}$. Due to Proposition 4, a unique chordal decomposition exists and we obtain

where

and we have used

The following theorem characterizes the properties of this optimization.

The theorem shows that, under suitable conditions, any local minimizer of (20) is globally optimal. This holds for any ES distribution with $\rho$ satisfying the assumptions. In particular, an important special case is chordal graphical models in MGGD:

Unfortunately, the recent result in [2] on $\rho(x)=\log(x)/2p$ does not satisfy our conditions and requires a different analysis. And there is no detailed proof to support the claim in [2] that this banded Tyler estimators converge to unique fixed points.

When the condition in Theorem 2 holds, the objective function is convex with respect to $\bm{C}$. Therefore, it can be numerically minimized via any general convex optimization solver. For example, in the MGGD case with $\beta\geq 0.5$ the problem can be expressed as

where $\bf{t}$ are auxiliary variables. This formulation with second order cone constraints can be easily solved using the popular CVX package [17].

Alternatively, the optimization can be addressed using a majorization - minimization (MM) technique, e.g. [5]. The method begins with an initial estimate $\hat{\bm{\Sigma}}_{0}\in\mathrm{S_{++}}(p)$ and updates it according the following iterations

where

$u(x)=\rho^{\prime}(x)$ and $C$ is chosen such that ${h(\bm{\Sigma}_{m},\bm{\Sigma}_{m})=L_{0}(\bm{\Sigma}_{m})}$. Since $\rho^{\prime\prime}(x)$ is continuous and nonnegative, we have

and $L_{0}(\bm{\Sigma}_{m})$ converges as $m\rightarrow\infty$. Each iteration step in (32)-(33) can be interpreted as a Gaussian graphical model optimization (the weights $u(\cdot)$ are constant with respect to the optimization variable). These minimizations have a simple closed form solution when $G(V,E)$ is chordal (see appendix). Thus, the proposed technique is very efficient for implementation in practice.

Finally, it is worth mentioning that our framework can also be extended to structure learning in MGGD graphical models. Structure learning, also known as covariance selection, considers the estimation of an inverse covariance which is known to be sparse but the sparsity pattern itself is unknown [13, 8, 29]. Adding a sparsity constraint usually destroys the convexity. Instead, the modern approach relies on a convex relaxation based on an L1 norm penalty. In our context, the problem is convex in the Cholesky factor rather than in the inverse covariance itself. This leads to the following convex minimization

where $\lambda$ is a regularization parameter, and $\|\bm{C}\|_{1}$ is a matrix version of the L1 norm, namely a sum over the absolute values of the elements in $\bm{C}$. It is important to emphasize that this approach is only applicable to chordal graphical models and it assumes that the perfect order of the variables is known a priori. Recent developments in high dimensional covariance estimation provide data-driven methods for identifying this order and structure [30, 28]. We leave this topic as a possible direction for future research.

Furthermore, our methodology can be easily extended to joint estimation of the covariance and the centering parameter $\mu$. For this purpose, consider the optimization of

where $\bm{x}=\bm{C}^{T}{\bm{\mu}}$ When the conditions in Theorem 2(a) hold, the objective function is jointly convex with respect to $(\bm{C},\bm{x})$, therefore any of its local minimizer is also its global minimizer. Its algorithm can be similarly addressed using a majorization - minimization (MM) technique as in (32)-(33).

In this section, we present results of numerical simulations. The purpose of these simulations is three fold: to validate our theoretical results using synthetic data, to provide insight on a few open theoretical questions using synthetic data and to demonstrate the usefulness of our theoretical models in a real world setting.

The first experiment considered inverse covariance estimation using synthetic data generated from a known MGGD distribution with $\beta\geq 0.5$ and a chordal structure. The theory in Section IV shows that, in this case, the MLE can be formulated as a convex optimization. To verify this result, we generated $n$ i.i.d. realizations of an MGGD with $p=10$, $\beta=0.5$ (corresponding to a classical multivariate Laplacian distribution), and a banded inverse covariance of width $b=4$. In particular, we used a Toeplitz inverse covariance with $1.0$ diagonal elements and $0.4$ off-diagonal elements within the main band. We estimated the unknown covariance using five estimators:

• •

  G: A Gaussian MLE with no prior knowledge of the structure, corresponding to the classical sample covariance.

• •

  BG: A Gaussian banded MLE as detailed in the appendix.

• •

  MGGD: An MGGD MLE with no prior knowledge of the structure via $30$ iterations of (14).

• •

  BMGGD1: An MGGD banded MLE via $30$ iterations of (32) and the subroutine in the appendix. This estimator is initialized with $\hat{\bm{\Sigma}}=\mathbf{I}$.

• •

  BMGGD2: An MGGD banded MLE via $30$ iterations of (32) and the subroutine in the appendix. This estimator is initialized with $\hat{\bm{\Sigma}}=\bm{\Sigma}$ (which is clearly impossible in practice).

In all estimators above we use $30$ as the number of iterations since in our simulations the fixed point algorithms (14) and (32) converge well after $30$ iterations.

Figure 1 shows the normalized mean squared Frobenius error in the covariance averaged over $10000$ independent simulations as a function of the number of samples $n$. It is easy to see the performance advantage of banded MGGD estimator. As expected, BMGGD1 and BMGGD2 converge to the identical fixed points irrespective of their initial condition.

The second experiment is very similar to the first except that $\beta=0.2$. Our theory assumes $\beta\geq 0.5$ and therefore does not hold for this value of shape parameter. Yet, according to [2] banded Tyler estimators do converge to unique fixed points, and these can be interpreted as the limit of MGGD when $\beta\rightarrow 0$. Our proposed fixed point iteration in (32) can be implemented with a small $\beta$. Therefore, it is interesting to examine its performance and we repeat the first experiment with $\beta=0.2$. The results are presented in Figure 2. Interestingly, BMGGD1 and BMGGD2 converged to identical fixed points irrespective of their initial condition. Thus, although we have no proof for this behavior, we conjecture that, in practice, the iteration can also be used for all MGGDs.

The third experiment focused on non-chordal structures. Our theory assumes a chordal sparsity pattern and it is interesting to see whether this assumption is indeed critical. Thus, we repeated the first experiment associated with MGGD $\beta=0.5$, but replaced the banded inverse covariance with a loopy structure. In particular, we constructed a two dimensional grid graph of size $p=3^{2}$ with the edges $(3,1)$, $(5,1)$, $(6,1)$, $(7,1)$, $(8,1)$, $(9,1)$, $(4,2)$, $(6,2)$, $(7,2)$, $(8,2)$, $(9,2)$, $(4,3)$, $(5,3)$, $(7,3)$, $(8,3)$, $(9,3)$, $(6,4)$, $(8,4)$, $(9,4)$, $(7,5)$, $(9,5)$, $(7,6)$, $(8,6)$, and $(9,7)$. All the non-zero off diagonal elements were assigned a constant value small enough to ensure that $\bm{\Sigma}^{-1}$ will be well conditioned. In contrast to the chordal case, MLE in general Gaussian Graphical models in does not satisfy a closed form solution. Instead, we solved (17) using the CVX optimization package [17]. The latter was used for implementing BG, and for the inner solution in each iteration of BMGGD1 and BMGGD2. The results averaged over $400$ simulations²²2Due to CVX these simulations are highly time consuming. are provided in Figure 3. Here too, the iterations converged to identical solutions and seem to be independent of their initial conditions.

The fourth experiment addressed the practical use of the BMGGD1 estimator in a real world example. Following [29] we considered the SONAR dataset from the UCI machine learning data repository. This dataset has 111 spectra from metal cylinders and 97 spectra from rocks, where each spectrum has 60 frequency band energy measurements. Quadratic discriminant test is used to classify the metals and the rocks. It requires the estimation of the covariance in both classes, and previous work demonstrated the advantage of BG over the classical sample covariance and the naive diagonal estimate. We repeated the experiment step by step and added BMGGD1. Specifically, we chose $\beta$ as the parameter within $\{0.5,0.6,\cdots,0.9,1.0\}$ that maximizes the MGGD likelihood, and used the same band which was used by BG (selected via 10 random splits with $1/3$ of the data for training and the validation likelihood). For the QDA test we applied the covariance of MGGD, which is $c(\beta)\bm{\Sigma}$, where $c(\beta)=2^{\frac{1}{\beta}}\frac{\Gamma(\frac{p+2}{2\beta})}{p\Gamma(\frac{p}{2\beta})}$. The test errors over a standard leave-one-out cross validation are provided in Table I. These results remained stable over different randomizations, and demonstrate the advantage of the proposed BMGGD1 framework.

In this paper, we consider covariance estimation in the multivariate generalized Gaussian distribution (MGGD). We proved that the MGGD negative log-likelihood is geodesically convex for $\beta>0$. In the sparsely constrained case, we proved that a simple change of variables can transform it into a convex function as long as $\beta\geq 0.5$ and the underlying graph is chordal. This means that any local solution of these minimization is globally optimal and the problems can be solved using standard descent methods. In practice, we observed this behavior also for smaller values of $\beta$. This agrees with a similar result on banded Tyler methods which can be interpreted as $\beta\rightarrow 0$. An interesting direction for future work is a rigorous analysis of these phenomenon when $0<\beta<0.5$. Another direction is relaxing the chordal assumption on the sparsity pattern. Here too, our numerical experience suggests that simple descent methods converge to the global solution and are independent of their initial conditions. Finally, as mentioned above, it in interesting to examine the problem of structure learning in MGGDs via sparsity enforcing priors.

In this appendix, we review a simple closed form solution for the MLE in chordal (decomposable) Gaussian graphical models [20, 34, 4]. The problem is

where

are the iteration weights which do not depend on $\bm{\Sigma}$. Using the chordal Cholesky decomposition $\bm{\Sigma}^{-1}=\bm{C}\bm{C}^{T}$, the problem is equivalent to

This minimization is completely separable and each column of $\bm{C}$ can be simply obtained by a linear regression. Let $\{\bm{Z}_{i}\}_{i=1}^{p}\subset\mathbb{R}^{n}$ be vectors consists of the $i$-th components of $\sqrt{\alpha_{1}}\bm{z}_{1},\sqrt{\alpha_{2}}\bm{z}_{2},\cdots,\sqrt{\alpha_{n}}\bm{z}_{n}$, $J_{i}=\{i+1\leq k\leq p:(i,k)\in E\}$, and assume that in the linear regression of $\bm{Z}_{i}$ with respect to $\{\bm{Z}_{j}\}_{j\in J_{i}}$, the parameter on $\bm{Z}_{j}$ is $\beta(i,j)$, and standard error is $r_{i}$. Then

and finally $\bm{C}=\bm{D}\bm{A}$.