Multilevel maximum likelihood estimation with application to covariance matrices

TL;DR

This paper demonstrates that restricting maximum likelihood estimation to a subspace can reduce asymptotic variance, improving covariance matrix fitting in high-dimensional data assimilation tasks.

Contribution

It introduces a multilevel maximum likelihood framework and applies it to spectral diagonal and sparse inverse covariance models, showing computational benefits.

Findings

Restricting to a subspace decreases asymptotic variance.

Smaller parameter sets reduce sampling noise.

Hierarchical approach improves covariance estimation in high dimensions.

Abstract

The asymptotic variance of the maximum likelihood estimate is proved to decrease when the maximization is restricted to a subspace that contains the true parameter value. Maximum likelihood estimation allows a systematic fitting of covariance models to the sample, which is important in data assimilation. The hierarchical maximum likelihood approach is applied to the spectral diagonal covariance model with different parameterizations of eigenvalue decay, and to the sparse inverse covariance model with specified parameter values on different sets of nonzero entries. It is shown computationally that using smaller sets of parameters can decrease the sampling noise in high dimension substantially.

Equations

$${\hat{\mathbf{\theta }}}_N=\arg \underset{\mathbf{\theta }}{\max }\mathcal{L}\left(\mathbf{\theta }|{\mathbb{X}}_N\right),\mathcal{L}\left(\mathbf{\theta }|{\mathbb{X}}_N\right)=\prod _{i=1}^N\mathcal{L}\left(\mathbf{\theta }|{\mathbf{X}}_i\right),\mathcal{L}\left(\mathbf{\theta }|\mathbf{x}\right)=f\left(\mathbf{x},\mathbf{\theta }\right),$$

$${\hat{\mathbf{\theta }}}_N=\arg \underset{\mathbf{\theta }}{\max }\ell \left(\mathbf{\theta }|{\mathbb{X}}_N\right),\ell \left(\mathbf{\theta }|{\mathbb{X}}_N\right)=\sum _{i=1}^N\ell \left(\mathbf{\theta }|{\mathbf{X}}_i\right),\ell \left(\mathbf{\theta }|\mathbf{x}\right)=\log f\left(\mathbf{x},\mathbf{\theta }\right).$$

$$\sqrt{N}({\hat{\mathbf{\theta }}}_N-{\mathbf{\theta }}^0)\stackrel{d}{\to }{\mathcal{N}}_p(0,Q_{{\mathbf{\theta }}^0}),\text{ as }N\to \infty ,$$

$$Q_{{\mathbf{\theta }}^0}=J_{{\mathbf{\theta }}^0}^{-1},J_{{\mathbf{\theta }}^0}=E\left({\nabla }_{\mathbf{\theta }}\ell ({\mathbf{\theta }}^0|\mathbf{X}{)}^{\top }{\nabla }_{\mathbf{\theta }}\ell ({\mathbf{\theta }}^0|\mathbf{X})\right),\mathbf{X}\sim f\left(\mathbf{x},{\mathbf{\theta }}^0\right).$$

$${\hat{\mathbf{\phi }}}_N=\arg \underset{\mathbf{\phi }}{\max }\ell \left(\mathbf{\phi }|{\mathbb{X}}_N\right),\ell \left(\mathbf{\phi }|{\mathbb{X}}_N\right)=\sum _{i=1}^N\ell \left(\mathbf{\phi }|{\mathbf{X}}_i\right),\ell \left(\mathbf{\phi }|\mathbf{x}\right)=\log f\left(\mathbf{x},\mathbf{\theta }\left(\mathbf{\phi }\right)\right),$$

$$\sqrt{N}({\hat{\mathbf{\phi }}}_N-{\mathbf{\phi }}^0)\stackrel{d}{\to }{\mathcal{N}}_k\left(0,Q_{{\mathbf{\phi }}^0}\right)\text{ as }N\to \infty ,$$

$$J_{{\mathbf{\phi }}^0}={\nabla }_{\mathbf{\phi }}\mathbf{\theta }({\mathbf{\phi }}^0{)}^{\top }{J}_{{\mathbf{\theta }}^0}{\nabla }_{\mathbf{\phi }}\mathbf{\theta }({\mathbf{\phi }}^0).$$

$${\nabla }_{\mathbf{\phi }}\ell (\mathbf{\phi }|\mathbf{X})={\nabla }_{\mathbf{\theta }}\ell (\mathbf{\theta }|\mathbf{X}){\nabla }_{\mathbf{\phi }}\mathbf{\theta }(\mathbf{\phi }),$$

$$J_{{\mathbf{\phi }}^0}$$

$$={\nabla }_{\mathbf{\phi }}\mathbf{\theta }({\mathbf{\phi }}^0{)}^{\top }E\left({\nabla }_{\mathbf{\theta }}\ell ({\mathbf{\theta }}^0|\mathbf{X}{)}^{\top }{\nabla }_{\mathbf{\theta }}\ell ({\mathbf{\theta }}^0|\mathbf{X})\right){\nabla }_{\mathbf{\phi }}\mathbf{\theta }({\mathbf{\phi }}^0)$$

$$={\nabla }_{\mathbf{\phi }}\mathbf{\theta }({\mathbf{\phi }}^0{)}^{\top }{J}_{{\mathbf{\theta }}^0}{\nabla }_{\mathbf{\phi }}\mathbf{\theta }({\mathbf{\phi }}^0).$$

$$\sqrt{N}(\mathbf{\theta }\left({\hat{\mathbf{\phi }}}_N\right)-{\mathbf{\theta }}^0)\stackrel{d}{\to }{\mathcal{N}}_p\left(0,Q_{\mathbf{\theta }\left({\mathbf{\phi }}^0\right)}\right)\text{ as }N\to \infty ,$$

$$Q_{\mathbf{\theta }\left({\mathbf{\phi }}^0\right)}={\nabla }_{\mathbf{\phi }}\mathbf{\theta }({\mathbf{\phi }}^0)J_{{\mathbf{\phi }}^0}^{-1}{\nabla }_{\mathbf{\phi }}\mathbf{\theta }({\mathbf{\phi }}^0{)}^{\top }={\nabla }_{\mathbf{\phi }}\mathbf{\theta }({\mathbf{\phi }}^0){\left({\nabla }_{\mathbf{\phi }}\mathbf{\theta }({\mathbf{\phi }}^0{)}^{\top }{J}_{{\mathbf{\theta }}^0}{\nabla }_{\mathbf{\phi }}\mathbf{\theta }({\mathbf{\phi }}^0)\right)}^{-1}{\nabla }_{\mathbf{\phi }}\mathbf{\theta }({\mathbf{\phi }}^0{)}^{\top }.$$

$$\sqrt{N}(\mathbf{\theta }({\hat{\mathbf{\psi }}}_N)-{\mathbf{\theta }}^0)\stackrel{d}{\to }{\mathcal{N}}_p\left(0,Q_{\mathbf{\theta }\left({\mathbf{\psi }}^0\right)}\right)\text{ as }N\to \infty ,$$

$$Q_{\mathbf{\theta }\left({\mathbf{\psi }}^0\right)}={\nabla }_{\mathbf{\psi }}\mathbf{\theta }({\mathbf{\psi }}^0)J_{{\mathbf{\psi }}^0}^{-1}{\nabla }_{\mathbf{\psi }}\mathbf{\theta }({\mathbf{\psi }}^0{)}^{\top }={\nabla }_{\mathbf{\psi }}\mathbf{\theta }({\mathbf{\psi }}^0){\left({\nabla }_{\mathbf{\psi }}\mathbf{\theta }({\mathbf{\psi }}^0{)}^{\top }{J}_{{\mathbf{\theta }}^0}{\nabla }_{\mathbf{\psi }}\mathbf{\theta }({\mathbf{\psi }}^0)\right)}^{-1}{\nabla }_{\mathbf{\psi }}\mathbf{\theta }({\mathbf{\psi }}^0{)}^{\top }.$$

$$Q_{\mathbf{\theta }\left({\mathbf{\phi }}^0\right)}\le Q_{\mathbf{\theta }\left({\mathbf{\psi }}^0\right)}.$$

$$E{\left|U\right|}^2=\frac{1}{N}\mathrm{Tr}{Q}_{\mathbf{\theta }\left({\mathbf{\phi }}^0\right)}\le \frac{1}{N}\mathrm{Tr}{Q}_{\mathbf{\theta }\left({\mathbf{\psi }}^0\right)}=E{\left|V\right|}^2,$$

$${\nabla }_{\mathbf{\phi }}\mathbf{\theta }\left({\mathbf{\phi }}^0\right)={\nabla }_{\mathbf{\psi }}\mathbf{\theta }\left({\mathbf{\psi }}^0\right){\nabla }_{\mathbf{\phi }}\mathbf{\psi }\left({\mathbf{\phi }}^0\right),$$

$$P_B$$

$$P_C$$

$$\mathrm{Range}{P}_B=\mathrm{Range}{A}^{1/2}B\subset \mathrm{Range}{A}^{1/2}C=\mathrm{Range}{P}_C.$$

$$\mathbf{X}\sim {\mathcal{N}}_n(0,D),$$

$$L\left(\Sigma |{\mathbb{X}}_N\right)=\frac{1}{{\left(\det \Sigma \right)}^{N/2}{\left(2\pi \right)}^{nN/2}}{e}^{-\frac{1}{2}\mathrm{Tr}\left({\Sigma }^{-1}{\mathbb{X}}_N{\mathbb{X}}_N^{\top }\right)}.$$

$${\hat{\Sigma }}_N=\frac{1}{N}\sum _{i=1}^N{\mathbf{X}}^{(i)}{\left({\mathbf{X}}^{(i)}\right)}^{\top }.$$

$$J^{(0)}(\mathrm{vec}(\Sigma ))=\frac{1}{2}{\Sigma }^{-1}\otimes {\Sigma }^{-1},$$

$$\tilde{\Sigma}_{N}=\left[\begin{array}[c]{cc}\tilde{\Sigma}_{11}&0\\ 0&0\end{array}\right],\quad\tilde{\Sigma}_{11}>0.$$

$$\underset{a\to 0^{+}}{\lim }L\left({\tilde{\Sigma }}_N+aI\right|{\mathbb{X}}_N)=\infty .$$

$${\hat{D}}_N^{\left(0\right)}=\mathrm{diag}\left({\hat{\Sigma }}_N\right).$$

$${\hat{\Sigma }}_N^u=\frac{1}{N-1}\sum _{i=1}^N{\mathbf{X}}^{(i)}{\left({\mathbf{X}}^{(i)}\right)}^{\top }$$

$${\left|{\hat{D}}_N^{\left(0\right)}-D\right|}_F\le {\left|{\hat{\Sigma }}_N^u-D\right|}_F$$

$$\ell (D|{\mathbb{X}}_N)=-\frac{N}{2}\log \left((2\pi {)}^n|D|\right)-\frac{1}{2}\sum _{k=1}^N{\left({\mathbf{X}}^{(k)}\right)}^{\top }{D}^{-1}{\mathbf{X}}^{(k)}$$

$${\hat{d}}_j=\frac{1}{N}\sum _{k=1}^N{\left(X_j^{(k)}\right)}^2,j=1,\dots ,n,$$

Full text

MULTILEVEL MAXIMUM LIKELIHOOD ESTIMATION WITH APPLICATION TO COVARIANCE MATRICES

Marie Turčičová*∗*

  

Jan Mandel*†*

  

Kryštof Eben*‡*

∗ Institute of Computer Science, Academy of Sciences of the Czech Republic

Pod Vodárenskou věží 271/2, 182 07 Praha 8, Czech Republic, and

Charles University in Prague, Faculty of Mathematics and Physics

Sokolovská 83, Prague 8, 186 75, Czech Republic

turcicova@cs.cas.cz † University of Colorado Denver, Denver, CO 80217-3364, USA, and

Institute of Computer Science, Academy of Sciences of the Czech Republic

Pod Vodárenskou věží 271/2, 182 07 Praha 8, Czech Republic

Jan.Mandel@ucdenver.edu ‡ Institute of Computer Science, Academy of Sciences of the Czech Republic

Pod Vodárenskou věží 271/2, 182 07 Praha 8, Czech Republic

eben@cs.cas.cz

Key Words: hierarchical maximum likelihood; nested parameter spaces; spectral diagonal covariance model; sparse inverse covariance model; Fisher information; high dimension.

ABSTRACT

The asymptotic variance of the maximum likelihood estimate is proved to decrease when the maximization is restricted to a subspace that contains the true parameter value. Maximum likelihood estimation allows a systematic fitting of covariance models to the sample, which is important in data assimilation. The hierarchical maximum likelihood approach is applied to the spectral diagonal covariance model with different parameterizations of eigenvalue decay, and to the sparse inverse covariance model with specified parameter values on different sets of nonzero entries. It is shown computationally that using smaller sets of parameters can decrease the sampling noise in high dimension substantially.

1 INTRODUCTION

Estimation of large covariance matrices from small samples is an important problem in many fields, including spatial statistics, genomics, and ensemble filtering. One of the prominent applications is data assimilation in meteorology and oceanography, where the dimension of state vector describing the atmosphere or ocean is in order of millions or larger. Every practically available sample is a small sample in this context, since a reasonable approximation of the full covariance can be obtained only with sample size of the order of the dimension of the problem (Vershynin, \APACyear2012). In practice, the sample covariance111In this paper, by sample covariance we mean the maximum likelihood estimate of covariance matrix using the norming constant $N$ as opposed to the unbiased estimate with norming constant $(N-1).$ is singular and polluted by spurious correlations. Nevertheless, it carries useful information (e.g. on covariances present in the actual atmospheric flow) and different techniques can be applied in order to improve the covariance model and its practical performance.

One common technique is shrinkage, that is, a linear combination of sample covariance and a positive definite target matrix, which prevents the covariance from being singular. The target matrix embodies some prior information about the covariance; it can be, e.g., unit diagonal or, more generally, positive diagonal (Ledoit \BBA Wolf, \APACyear2004). See, e.g., Schäfer \BBA Strimmer (\APACyear2005) for a survey of such shrinkage approaches. Shrinkage of sample covariance towards a fixed covariance matrix based on a specific model and estimated from historical data (called background covariance) was used successfully in meteorology (Hamill \BBA Snyder, \APACyear2000; Wang \BOthers., \APACyear2008). This approach is justified as one which combines actual (called flow-dependent) and long-term average (called climatologic) information on spatial covariances present in the 3D meteorological fields.

Another approach to improving on the sample covariance matrix is localization by suppressing long-range spurious correlations, which is commonly done by multiplying the sample covariance matrix term by term by a gradual cutoff matrix (Buehner \BBA Charron, \APACyear2007; Furrer \BBA Bengtsson, \APACyear2007) to suppress off-diagonal entries. The extreme case, when only the diagonal is left, is particularly advantageous in the spectral domain, as the covariance of a random field in Fourier space is diagonal if and only if the random field in cartesian geometry is second order stationary, i.e., the covariance between the values at two points depends only on their distance vector. Alternatively, diagonal covariance in a wavelet basis provides spatial variability as well (Pannekoucke \BOthers., \APACyear2007). Spectral diagonal covariance models were successfully used in operational statistical interpolation in meteorology in spherical geometry (Parrish \BBA Derber, \APACyear1992), and versions of Ensemble Kalman Filter (EnKF) were developed which construct diagonal covariance in Fourier or wavelet space in every update step of the filter at low cost, and can operate successfully with small ensembles (Mandel \BOthers., \APACyear2010; Beezley \BOthers., \APACyear2011; Kasanický \BOthers., \APACyear2015).

Sparse covariance models, such as the spectral diagonal, allow a compromise between realistic assumptions and cheap computations. Another covariance model taking advantage of sparsity is a Gauss-Markov Random Field (GMRF), based on the fact that conditional independence of variables implies zero corresponding elements in the inverse of the covariance matrix (Rue \BBA Held, \APACyear2005), which leads to modeling the covariance as the inverse of a sparse matrix.

However, both spectral diagonal and sparse inverse covariance models have a large number of parameters, namely all terms of the sparse matrix (up to symmetry) which are allowed to attain nonzero values. This results in overfitting and significant sampling noise for small samples. Therefore, it is of interest to reduce the number of parameters by adopting additional, problem-dependent assumptions on the true parameter values.

The principal result of this paper is the observation that if parameters are fitted as the Maximum Likelihood Estimator (MLE) and the additional assumptions are satisfied by the true parameters, then the estimate using fewer parameters is asymptotically more accurate, and often very significantly so even for small samples..

The paper is organized as follows. In Sec. 2, we provide a brief statement of MLE and its asymptotic variance. In Sec. 3, we use the theory of maximum likelihood estimation to prove that for any two nested subspaces of the parametric space containing the true parameter, the asymptotic covariance matrix of the MLE is smaller for the smaller parameter space. These results hold for a general parameter and, in the special case of MLE for covariance matrices we do not need any invertibility assumption. The applications to estimation of covariance matrices by spectral diagonal and GMRF are presented in Sec. 4, and Sec. 5 contains computational illustrations. A comparison of the performance of MLE for parametric models and of related shrinkage estimators is in Sec. 6.

2 ASYMPTOTIC VARIANCE OF THE MAXIMUM LIKELIHOOD ESTIMATOR

First, we briefly review some standard results for reference. Suppose $\mathbb{X}_{N}=\left[\bm{X}_{1},\ldots,\bm{X}_{N}\right]$ is a random sample from a distribution on $\mathbb{R}^{n}$ with density $f\left(\bm{x},\bm{\theta}\right)$ with unknown parameter vector $\bm{\theta}$ in a parameter space $\Theta\subset\mathbb{R}^{p}$. The maximum likelihood estimate $\bm{\hat{\theta}}_{N}$ of the true parameter $\bm{\theta}^{0}$ is defined by maximizing the likelihood

[TABLE]

or, equivalently, maximizing the log likelihood

[TABLE]

We adopt the usual assumptions that (i) the true parameter $\bm{\theta}^{0}$ lies in the interior of $\Theta$, (ii) the density $f$ determines the parameter $\bm{\theta}$ uniquely in the sense that $f(\bm{x},\bm{\theta}_{1})=f(\bm{x},\bm{\theta}_{2})$ a.s. if and only if $\bm{\theta}_{1}=\bm{\theta}_{2}$, and (iii) $f\left(\bm{x},\bm{\theta}\right)$ is a sufficiently smooth function of $\bm{x}$ and $\bm{\theta.}$ Then the error of the estimate is asymptotically normal,

[TABLE]

where

[TABLE]

The matrix $J_{\bm{\theta}^{0}}$ is called the Fisher information matrix for the parameterization $\bm{\theta}^{0}$. Here, $\bm{X}$, $\bm{x}$, and $\bm{\theta}$ are columns, while the gradient $\nabla_{\bm{\theta}}\ell$ of $\ell$ with respect to the parameter $\bm{\theta}$ is a row vector, which is compatible with the dimensioning of Jacobi matrices below. The mean value in (3) is taken with respect to $\bm{X}$, which is the only random quantity in (3). Cf., e.g., (Lehmann \BBA Casella, \APACyear1998, Theorem 5.1) for details.

3 NESTED MAXIMUM LIKELIHOOD ESTIMATORS

Now suppose that we have an additional information that the true parameter $\bm{\theta}^{0}$ lies in a subspace of $\Theta$, which is parameterized by $k\leq p$ parameters $(\varphi_{1}$,…, $\varphi_{k})^{\top}=\bm{\varphi}$. Denote by $\nabla_{\bm{\varphi}}\bm{\theta}(\bm{\varphi})$ the $p\times k$ Jacobi matrix with entries $\frac{\partial\theta_{i}}{\partial\varphi_{j}}$. In the next theorem, we derive the asymptotic covariance of the maximum likelihood estimator for $\bm{\varphi}$,

[TABLE]

based on the asymptotic covariance of $\bm{\theta}$ in (2).

Theorem 1

Assume that the map $\bm{\varphi}\mapsto\bm{\theta}(\bm{\varphi})$ is one-to-one from $\Phi\subset\mathbb{R}^{k}$ to $\Theta$, the map $\bm{\varphi}\mapsto\bm{\theta}(\bm{\varphi})$ is continuously differentiable, $\nabla_{\bm{\varphi}}\bm{\theta}(\bm{\varphi})$ is full rank for all $\bm{\varphi}\in\Phi$, and $\bm{\theta}^{0}=\bm{\theta}(\bm{\varphi}^{0})$ with $\bm{\varphi}^{0}$ in the interior of $\Phi$. Then,

[TABLE]

where $Q_{\bm{\varphi}^{0}}=J_{\bm{\varphi}^{0}}^{-1}$, with $J_{\bm{\varphi}^{0}}$ the Fisher information matrix of the parameterization $\bm{\varphi}$ given by

[TABLE]

Proof. From (3) and the chain rule

[TABLE]

we have

[TABLE]

The asymptotic distribution (5) is now (2) applied to $\bm{\varphi}$.

When the parameter $\bm{\theta}$ is the quantity of interest in an application, it is useful to express the estimate and its variance in terms of the original parameter $\bm{\theta}$ rather than the subspace parameter $\bm{\varphi}$.

Corollary 2

Under the assumptions of Theorem 1,

[TABLE]

where

[TABLE]

Proof. The lemma follows from (5) by the delta method (Rao, \APACyear1973, p. 387), since the map $\bm{\varphi\mapsto\theta}(\bm{\varphi})$ is continuously differentiable.

Remark 3

The matrix $Q_{\bm{\theta}\left(\bm{\varphi}^{0}\right)}$ is singular, so it cannot be written as the inverse of another matrix, but it can be understood as the inverse $J_{\bm{\theta}(\bm{\varphi}^{0})}^{-1}$ of the Fisher information matrix for $\bm{\varphi}$, embedded in the larger parameter space $\Theta$.

Suppose that $\bm{\psi}$ is another parameterization which satisfies the same assumption as $\bm{\varphi}$ in Theorem 1: the map $\bm{\psi}\mapsto\bm{\theta}(\bm{\psi})$ is one-to-one from $\Psi\subset\mathbb{R}^{m}$, $k\leq m\leq p,$ to $\Theta$, $\bm{\psi}\mapsto\bm{\theta}(\bm{\psi})$ is continuously differentiable, $\nabla_{\bm{\psi}}\bm{\theta}(\bm{\psi})$ is full rank for all $\bm{\psi}\in\Psi$, and $\bm{\theta}^{0}=\bm{\theta}(\bm{\psi}^{0})$, where $\bm{\psi}^{0}$ is in the interior of $\Psi$. Then, similarly as in (7), we have also

[TABLE]

where, as in (8),

[TABLE]

The next theorem shows that when we have two parameterizations $\bm{\varphi}$ and $\bm{\psi}$ which are nested, then the smaller parameterization has smaller or equal asymptotic covariance than the larger one. For symmetric matrices $A$ and $B$, $A\leq B$ means that $A-B$ is positive semidefinite.

Theorem 4

Suppose that $\bm{\varphi}$ and $\bm{\psi}$ satisfy the assumptions in Theorem 1, and there exists a differentiable mapping $\bm{\varphi\mapsto\psi}$ from $\Phi$ to $\Psi$, such that $\bm{\varphi}^{0}\bm{\mapsto\psi}^{0}$. Then,

[TABLE]

In addition, if $U\sim\mathcal{N}_{p}\left(\bm{0},Q_{\bm{\theta}\left(\bm{\varphi}^{0}\right)}\right)$ and $V\sim\mathcal{N}_{p}\left(\bm{0},Q_{\bm{\theta}\left(\bm{\psi}^{0}\right)}\right)$ are random vectors with the asymptotic distributions of the estimates $\bm{\theta}\left(\hat{\bm{\varphi}}_{N}\right)$ and $\bm{\theta}\left(\hat{\bm{\psi}}_{N}\right)$, then

[TABLE]

where $\left|V\right|=\left(V^{\top}V\right)^{1/2}$ is the standard Euclidean norm in $\mathbb{R}^{p}$.

Proof. Denote $A=J_{\bm{\theta}^{0}}$, $B=\nabla_{\bm{\varphi}}\bm{\theta}(\bm{\varphi}^{0})$, $C=\nabla_{\bm{\psi}}\bm{\theta}(\bm{\psi}^{0})$. From the chain rule,

[TABLE]

we have that $B=C\nabla_{\bm{\varphi}}\bm{\psi}\left(\bm{\varphi}^{0}\right)$, and, consequently, $\operatorname*{Range}B\subset\operatorname*{Range}C$. Define

[TABLE]

The matrices $P_{B}$ and $P_{C}$ are symmetric and idempotent, hence they are orthogonal projections. In addition,

[TABLE]

Consequently, $P_{B}\leq P_{C}$ holds from standard properties of orthogonal projections, and (11) follows.

To prove (12), note that for random vector $X$ with $\operatorname{E}X=0$ and finite second moment, $\operatorname{E}\left|X\right|^{2}=\operatorname*{Tr}\operatorname{Cov}X$ from Karhunen-Loève decomposition and Parseval identity. The proof is concluded by using the fact that for symmetric matrices, $A\leq B$ implies $\operatorname*{Tr}A\leq\operatorname*{Tr}B$, cf. e.g., Carlen (\APACyear2010).

Remark 5

In the practically interesting cases when there is a large difference in the dimensions of the parameters $\bm{\varphi}$ and $\bm{\psi}$, many eigenvalues in the covariance of the estimation error become zero. The computational tests in Sec. 5 show that the resulting decrease of the estimation error can be significant.

4 APPLICATION: NESTED COVARIANCE MODELS

Models of covariance (e.g., of the state vector in a numerical weather prediction model) and the quality of the estimated covariance are one of the key components of data assimilation algorithms. High dimension of the problem usually prohibits working with the covariance matrix explicitly. In ensemble filtering methods, this difficulty may be circumvented by working directly with the original small sample like in the classical Ensemble Kalman filter. This, however, effectively means using the sample covariance matrix with its rank deficiency and spurious correlations. Current filtering methods use shrinkage and localization as noted above, and ad hoc techniques for dimension reduction.

A reliable way towards effective filtering methods lies in introducing sparsity into covariance matrices or their inverses by means of suitable covariance models. The results of previous section suggest that it is beneficial to choose parsimonious models, and indeed, in practical application we often encounter models with a surprisingly low number of parameters.

A large class of covariance models which encompass sparsity in an efficient manner arises from Graphical models (Lauritzen, \APACyear1996) and Gaussian Markov Random Fields (GMRF), (Rue \BBA Held, \APACyear2005), where a special structure of inverse covariance is assumed. In the area of GMRF, nested covariance models arise naturally. If, for instance, we consider a GMRF on a rectangular mesh, each gridpoint may have 4, 8, 12, 20 etc. neighbouring points which have nonzero corresponding element in the inverse covariance matrix. Thus, a block band-diagonal structure in the inverse covariance arises (Ueno \BBA Tsuchiya, \APACyear2009). The results of Section 3 apply for this case and we shall illustrate them in the simulation study of Section 5.

Finally, variational assimilation methods, which dominate today’s practice of meteorological services, usually employ a covariance model based on a series of transformations leading to independence of variables (Bannister, \APACyear2008; Michel \BBA Auligné, \APACyear2010). At the end, this results in an estimation problem for normal distribution with a diagonal covariance matrix.

For both ensemble and variational methods, any additional knowledge can be used to improve the estimate of covariance. Second-order stationarity leads to diagonality in spectral space, diagonality in wavelet space is often a legitimate assumption (e.g., Pannekoucke \BOthers., \APACyear2007) and we shall treat the diagonal case in more detail. Suppose

[TABLE]

where $\bm{X}$ denotes the random field after the appropriate transform and $D$ is a diagonal matrix. It is clear that estimating $D$ by the full sample covariance matrix (what would be the case when using the classical EnKF) is ineffective in this situation and it is natural to use only the diagonal part of the sample covariance. In practice, the resulting diagonal matrix may still turn out to be noisy (Kasanický \BOthers., \APACyear2015), and further assumptions like a certain type of decay of the diagonal entries may be realistic.

In what follows we briefly introduce the particular covariance structures, state some known facts on full and diagonal covariance, propose parametric models for the diagonal and compute corresponding MLE.

4.1 Sample covariance

The top-level parameter space $\Theta$ consists of all symmetric positive definite matrices, resulting in the parameterization $\Sigma$ with $\frac{n\left(n+1\right)}{2}$ independent parameters. The likelihood of a sample $\mathbb{X}_{N}=\left[\bm{X}^{(1)},\ldots,\bm{X}^{(N)}\right]$ from $\mathcal{N}_{n}(\bm{0},\Sigma)$ is

[TABLE]

If $N\geq n$, it is well known (e.g. Muirhead (\APACyear2005), p. 83) that the likelihood is maximized at what we call here sample covariance matrix

[TABLE]

The Fisher information matrix of the sample covariance estimator is (Magnus \BBA Neudecker, \APACyear2007, p. 356)

[TABLE]

where $\otimes$ stands for the Kronecker product and $\operatorname{vec}$ is an operator that transforms a matrix into a vector by stacking the columns of the matrix one underneath the other. This matrix has dimension $n^{2}\times n^{2}$.

Remark 6

If $\hat{\Sigma}_{N}$ is singular, $L\left(\hat{\Sigma}_{N}|\mathbb{X}_{N}\right)$ cannot be evaluated because that requires the inverse of $\hat{\Sigma}$. Also, in this case the likelihood $L\left(\Sigma|\mathbb{X}_{N}\right)$ is not bounded above on the set of all $\Sigma>0$, thus the maximum of $L\left(\Sigma|\mathbb{X}_{N}\right)$ does not exist on that space. To show that, consider an orthonormal change of basis so that the vectors in $\operatorname*{span}\left(\mathbb{X}_{N}\right)$ come first, write vectors and matrices in the corresponding $2\times 2$ block form, and let

[TABLE]

Then $\lim_{a\rightarrow 0^{+}}\mathbb{X}_{N}^{\top}\left(\tilde{\Sigma}_{N}+aI\right)^{-1}\mathbb{X}_{N}$ exists, but $\lim_{a\rightarrow 0^{+}}\det\left(\tilde{\Sigma}_{N}+aI\right)=0$, thus

[TABLE]

Note that when the likelihood is redefined in terms of the subspace $\operatorname*{span}\left(\mathbb{X}_{N}\right)$ only, the sample covariance can be obtained by maximization on the subspace (Rao, \APACyear1973, p. 527).

When the true covariance is diagonal ($\Sigma\equiv D$, cf. (13)), a significant improvement can be achieved by setting the off-diagonal terms of sample covariance to zero,

[TABLE]

It is known that using only the diagonal of the unbiased sample covariance

[TABLE]

results in smaller (or equal) Frobenius norm of the error pointwise,

[TABLE]

cf. Furrer \BBA Bengtsson (\APACyear2007) for the case when the mean is assumed to be known like here, and Kasanický \BOthers. (\APACyear2015) for the unbiased sample covariance and unknown mean.

4.2 Diagonal covariance

The parameter space $\Theta_{1}$ consisting of all diagonal matrices with positive diagonal, with $n$ parameters $\bm{d}=\left(d_{1},\ldots,d_{n}\right)^{\top}$, can be viewed as a simple class of models for either covariance or its inverse. The log-likelihood function for $D=\operatorname{diag}(d_{1},\ldots,d_{n})$ with a given random sample $\mathbb{X}_{N}=\left[\bm{X}^{(1)},\ldots,\bm{X}^{(N)}\right]$ from $\mathcal{N}_{n}\left(\bm{0},D\right)$ is

[TABLE]

and has its maximum at

[TABLE]

where $X_{j}^{(k)}$ denotes the $j$-th entry of $\bm{X}^{(k)}$. The sum of squares $S_{j}^{2}=\sum_{k=1}^{N}\left(X_{j}^{(k)}\right)^{2}$ is a sufficient statistic for the variance $d_{j}$. Thus, we get the maximum likelihood estimator

[TABLE]

It is easy to compute the Fisher information matrix explicitly,

[TABLE]

which is an $n\times n$ matrix and gives the asymptotic covariance of the estimation error

[TABLE]

from (2).

4.3 Diagonal covariance with prescribed decay by 3 parameters

A more specific situation appears when we have an additional information that the matrix $D$ is not only diagonal, but its diagonal entries have a prescribed decay. For instance, this decay can be governed by a model of the form $d_{i}=((c_{1}+c_{2}h_{i})f_{i}(\alpha))^{-1}$, $i=1,\ldots,n$, where $c_{1},c_{2}$ and $\alpha$ are unknown parameters, $h_{1},\ldots,h_{n}$ are known positive numbers, and $f_{1},\ldots,f_{n}$ are known differentiable functions. For easier computation it is useful to work with $\tau_{i}=\frac{1}{d_{i}}=(c_{1}+c_{2}h_{i})f_{i}(\alpha)$. Maximum likelihood estimators for $c_{1},c_{2}$, and $\alpha$ can be computed effectively from the likelihood

[TABLE]

by using the chain rule. It holds that

[TABLE]

Setting this derivative equal to zero we get

[TABLE]

Analogously,

[TABLE]

so the equation for estimating the parameter $c_{2}$ is

[TABLE]

Similarly,

[TABLE]

and setting the derivative to zero, we get

[TABLE]

The maximum likelihood estimator for $D$ is then given by

[TABLE]

where $(\hat{c}_{1},\hat{c}_{2},\hat{\alpha})$ is the solution of the system (21), (22), (23). This expression corresponds to searching a maximum likelihood estimator of $D$ in the subspace $\Theta_{3}\subset\Theta_{1}\subset\Theta$ formed by diagonal matrices $\operatorname{diag}\{((c_{1}+c_{2}h_{i})f_{i}(\alpha))^{-1},i=1,\ldots,n\}$.

For completeness, the asymptotic covariance of the estimation error about

[TABLE]

contained in $\mathbb{X}_{N}$ is

[TABLE]

from (8), where the Fisher information matrix $J_{c_{1},c_{2},\alpha}$ is the $3\times 3$ matrix

[TABLE]

and

[TABLE]

4.4 Diagonal covariance with prescribed decay by 2 parameters

We may consider a more specific model for diagonal elements with two parameters: $d_{i}=(cf_{i}(\alpha))^{-1}$, i.e. $\tau_{i}=cf_{i}(\alpha)$, $i=1,\ldots,n$, where $c$ and $\alpha$ are unknown parameters. Maximum likelihood estimators for $c$ and $\alpha$ can be computed similarly as in the previous case. The estimating equations have the form

[TABLE]

which can be rearranged to

[TABLE]

Equation (27) is an implicit formula for estimating $\alpha$. Its result can be used for estimating $c$ through (26). The maximum likelihood estimator for $D$ is then given by

[TABLE]

where $\hat{c}$ and $\hat{\alpha}$ are MLEs of $c$ and $\alpha$. It corresponds to searching a maximum likelihood estimator of $D$ in the subspace $\Theta_{2}\subset\Theta_{3}\subset\Theta_{1}\subset\Theta$ formed by diagonal matrices $\operatorname{diag}\left\{(cf_{i}(\alpha))^{-1},i=1,\ldots,n\right\}$. Of course, the estimator $\hat{D}^{(2)}$ does not have “larger” variance than $\hat{D}^{(3)}$.

The covariance of the asymptotic distribution of the parameters $d_{1},\ldots,d_{n}$ is

[TABLE]

from (8), where Fisher information matrix at $D=\operatorname{diag}\{d_{i}(c,\alpha),i=1,\ldots,n\}$ is the $2\times 2$ matrix

[TABLE]

and $\bm{d}(c,\alpha)=\left[d_{1}(c,\alpha),\ldots,d_{n}(c,\alpha)\right]^{\top}=\left[(cf_{1}(\alpha))^{-1},\ldots,(cf_{n}(\alpha))^{-1}\right]^{\top}.$

4.5 Sparse inverse covariance and GMRF

In the GMRF method for fields on a rectangular mesh, we assume that a variable on a gridpoint is conditionally independent on the rest of the gridpoints, given values on neighboring gridpoints. It follows that nonzero entries in the inverse of the covariance matrix can be only between neighbor gridpoints. We start with 4 neighbors (up, down, right, left), and adding neighbors gives rise to a sequence of nested covariance models. If the columns of the mesh are stacked vertically, their inverse covariance matrix will have a band-diagonal structure.

The inverse covariance model fitted by MLE was introduced by Ueno \BBA Tsuchiya (\APACyear2009) and applied on data from oceanography. The corresponding Fisher information matrix may be found as the negative of the Hessian matrix (Ueno \BBA Tsuchiya, \APACyear2009, eq. (C17)).

5 COMPUTATIONAL STUDY

In Section 3, we have shown that in the sense of asymptotic variance and second moment (mean-squared) error, the maximum likelihood estimator computed in a smaller space containing the true parameter is more (or equally) precise. For small samples, we illustrate this behavior by means of simulations.

5.1 Simulation of simple GMRF

We first show that in the case of GMRF with four neighbors per gridpoint, adding dependencies (parameters) which are not present brings a loss of precision of the MLE. Using the sample covariance in this case causes a substantial error.

We have generated an ensemble of realizations of a GMRF with dimensions $10\times 10$ (resulting in $n=100$) and inverse covariance structure as in Fig. 1. The values on the diagonals of the covariance matrix have been set to constant, since we assume the correlation with left and right neighbor to be identical, as well as the correlation with upper and lower neighbor (by symmetry of the covariance matrix and isotropy in both directions of the field, but different correlation in each direction). This leads to a model with 3 parameters for 4 neighbors, 5 parameters for 8 neighbors and 7 parameters for 12 neighbors,

The covariance structure of $\Sigma^{-1}$ with 4 neighbors was set as “truth” and random samples were generated from $\mathcal{N}_{n}(\bm{0},\Sigma)$ with sample sizes $N=10,15,20,\ldots,55$. The values on first, second and tenth diagonal have been set as 5, -0.2 and 0.5. For each sample, we computed successively the MLE with 3, 5 and 7 unknown parameters numerically by Newton’s method, as described in Ueno \BBA Tsuchiya (\APACyear2009).

The difference of each estimator from the true matrix $\Sigma$ was measured in the Frobenius norm, which is the same as the Euclidean norm of a matrix written as one long vector. In order to reduce the sampling error, 50 simulations of the same size were generated and the mean of squared Frobenius norm was computed. The results can be found in Fig. 2.

As expected, the MLE with 3 parameters outperforms the estimates with 5 and 7 parameters and the Frobenius norm for sample covariance stays one order worse than all parametric estimates.

5.2 Simulation of fields with diagonal covariance

The simulation for spectral diagonal covariance was carried out in a similar way. First, a diagonal matrix $D$ was prepared, whose diagonal entries decay according to the model $d_{i}=\frac{1}{c}e^{\alpha\lambda_{i}},i=1,\ldots,n$, where $c$ and $\alpha$ are parameters and $\lambda_{i}$ are the eigenvalues of Laplace operator in two dimensions on $10\times 10$ nodes (so again $n=100$). Such models are useful in modeling smooth random fields, e.g., in meteorology. Then, random samples were generated from $\mathcal{N}_{n}(\bm{0},D)$ with sample sizes $N=5,\ldots,20$. For each sample, four covariance matrix estimators were computed:

• •

sample covariance matrix $\hat{\Sigma}_{N}$, cf. (14)

• •

diagonal part $\hat{D}^{\left(0\right)}$ of the sample covariance matrix, cf. (15)

• •

MLE $\hat{D}^{\left(1\right)}$ in the space of diagonal matrices, cf. (18)

• •

MLE $\hat{D}^{\left(3\right)}=\operatorname{diag}\{(\hat{c}_{1}-\hat{c}_{2}\lambda_{i})^{-1}e^{\hat{\alpha}\lambda_{i}},i=1,\ldots,n\}$ with 3 parameters $c_{1},c_{2}$ and $\alpha$, cf. (24).

• •

MLE $\hat{D}^{\left(2\right)}=\operatorname{diag}\{\hat{c}^{-1}e^{\hat{\alpha}\lambda_{i}},i=1,\ldots,n\}$ with 2 parameters $c$ and $\alpha$, cf. (28).

Let us briefly discuss the choice of the covariance model $d_{i}=\frac{1}{c}e^{\alpha\lambda_{i}}$. We decided to carry out the simulation with a second-order stationary random field, whose covariance can be diagonalized by the Fourier transform. This transform is formed by the eigenvectors of the Laplace operator. Hence, it is reasonable to model the diagonal terms of this covariance matrix (i.e. the covariance eigenvalues) by some function of eigenvalues of the Laplace operator. This function needs to have a sufficiently fast decay in order to fulfil the necessary condition for the proper covariance (the so-called trace class property, e.g., Kuo (\APACyear1975)). Exponential decay is used, e.g., in Mirouze \BBA Weaver (\APACyear2010). Another possible choice of a covariance model is a power model, where the eigenvalues of the covariance are assumed to be a negative power of $-\lambda_{i},i=1,\ldots,n$, e.g., Berner \BOthers. (\APACyear2009); Gaspari \BOthers. (\APACyear2006); Simpson \BOthers. (\APACyear2012).

The difference of each estimator from the true matrix $D$ was measured in the Frobenius norm again. To reduce the sampling noise, 50 replications have been done for each sample size and the mean of squared Frobenius norm can be found in Fig. 3.

For the diagonal MLE, given by (18), (24), and (28), we can expect from (12) that these estimators should satisfy asymptotically

[TABLE]

even if convergence in distribution does not imply convergence of moments without additional assumptions. This conjecture can be supported by a comparison of Figures 5 and 4, where we observe the same decay. From the nesting, we know that

[TABLE]

and we can expect that the Frobenius norm should decrease for more restrictive models, that is,

[TABLE]

which is confirmed by the simulations (see Figure 4, resp. 5).

The comparisons (32) of the Frobenius norm of the error in the mean squared complement the pointwise comparison (16) between the sample covariance and its diagonal. Relying on MLE for that comparison is not practical, because the sample size of interest here is $N<n$, and, consequently, $\hat{\Sigma}_{N}$ is singular and cannot be cast as MLE with an accompanying Fisher information matrix, cf. Remark 6. But it is evident that for small sample sizes, estimators computed in the proper subspace perform better. Hence, the hierarchical order seems to hold even when $N<n$.

6 COMPARISON WITH REGULARIZATION METHODS

In the previous sections, we pointed out the advantages of using low-parametric models for estimating a covariance matrix using a small sample. As mentioned in the Introduction, there is another large class of estimating methods for high-dimensional covariance matrices: shrinkage estimators. The principle of these methods is to move the sample covariance towards a target matrix that possesses some desired properties (e.g., full rank, proper structure). This can be seen as a convex combination of the sample covariance matrix $\hat{\Sigma}_{N}$ and the so called target matrix $T$:

[TABLE]

One of the simplest shrinkage estimators has the form of (33) with the target matrix equal to identity, which results in shrinking all sample eigenvalues with the same intensity towards their mean value. Ledoit \BBA Wolf (\APACyear2004) derived the optimal shrinkage parameter $\gamma$ to minimize the squared Frobenius loss

[TABLE]

The comparison of this estimator with the maximum likelihood estimator $\hat{D}^{(2)}$ was accomplished by a simulation with identical setting as in Section 5. The results are shown in Fig. 6. For reference, the sample covariance $\hat{\Sigma}_{N}$ and its diagonal $\hat{D}^{(0)}$ are also added.

Another regularization method is described in Won \BOthers. (\APACyear2013). They consider a type of covariance estimator, where the regularization effect is achieved by bounding the condition number of the estimate by a regularization parameter $\kappa_{max}$. Since the condition number is defined as a ratio of the largest and smallest eigenvalue, this method corrects for overestimation of the largest eigenvalues and underestimation of the small eigenvalues simultaneously. The resulting estimator is called a condition-number-regularized covariance estimator and it is formulated as the maximum likelihood estimator restricted on the subspace of matrices with condition number bounded by $\kappa_{max}$, i.e.

[TABLE]

where $\lambda_{max}(\Sigma)$, resp. $\lambda_{min}(\Sigma)$, is the largest, resp. the smallest, eigenvalue of the covariance matrix $\Sigma$. An optimal $\kappa_{max}$ is selected by maximization of the expected likelihood, which is approximated by using $K$-fold cross-validation. The authors proved that $\kappa_{max}$ selected in this way is a consistent estimator for the true condition number (i.e. the condition number of $D$). Therefore, the idea of this method is to search a MLE in a subspace defined by covariance matrices with condition number smaller or equal to the true condition number. The form of the resulting covariance estimator together with the details of the computational process is provided in Won \BOthers. (\APACyear2013). In Fig. 6, we can see the performance of this estimator (denoted as cond-num-regularization) in comparison of other methods.

The shrinkage estimator $\hat{\Sigma}_{S}$ and the condition-number-regularized estimator result in non-diagonal matrices, which in our case predetermines them to perform worse than the diagonal estimator $\hat{D}^{(0)}$. However, we have to note that performance of these methods strongly depends on the particular form of the true covariance matrix $D$. In the case when the decrease of the true eigenvalues is less rapid, both methods may provide better results than the diagonal of sample covariance. The performance of $\hat{\Sigma}_{S}$ could be possibly improved by choosing a different target matrix that is closer to reality but such a study is out of the scope of this paper.

It is seen from Fig. 6 that the condition-number-regularized estimator provides more precise estimates than the sample covariance $\hat{\Sigma}_{N}$, as expected. This is in accordance with the preceding theory and empirical findings about the higher precision of estimators from a smaller parametric subspace (the corresponding parametric subspace consists of matrices with the condition number smaller or equal to $\kappa_{max}$). If, however, the theoretical condition number is very large as in our case, the method has a problem in estimating this number and its performance is limited.

Both regularization estimators perform well against sample covariance, but the setting of our simulation is less favourable for them. Neither of them can compete with the maximum likelihood estimator found in the true small subspace of diagonal matrices with proper decay.

7 CONCLUSIONS

Our main aim was to point out the significant advantage resulting from computing the MLE of the covariance matrix in a proper parameter subspace, especially in the high-dimensional setting, when the available sample has small size relative to the dimension of the problem. This subspace can be formed, e.g., by a parametric model for covariance eigenvalues or for a diagonal matrix resulting from a suitable set of transformations.

We provided theoretical results on asymptotic comparison of covariance matrices of each estimator for multivariate normal distribution, where we can lean on the well-developed maximum likelihood theory. The situation for small samples was illustrated by means of a simulation. We consider two-parametric models for the covariance eigenvalues based on the eigenvalues of Laplace operator. In practice, the proper model/subspace can be inferred from historical data.

Using a properly specified model, one can reach a significant improvement in performance, which can have a positive impact on the subsequent tasks like data assimilation and prediction.

ACKNOWLEDGEMENTS

This work was partially supported by the the Czech Science Foundation (GACR) under grant 13-34856S and by the U.S. National Science Foundation under grants DMS-1216481 and ICER-1664175.