Multilevel maximum likelihood estimation with application to covariance matrices
Marie Turčičová Jan Mandel Kryštof Eben





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 RepublicPod Vodárenskou věží 271/2, 182 07
Praha 8, Czech Republic, and Charles University in Prague, Faculty of
Mathematics and PhysicsSokolovská 83, Prague 8, 186 75, Czech
Republicturcicova@cs.cas.cz 3mm ^†
University of Colorado Denver, Denver, CO 80217-3364, USA, and Institute of Computer Science, Academy of Sciences of the Czech
RepublicPod Vodárenskou věží 271/2, 182 07 Praha 8,
Czech Republic Jan.Mandel@ucdenver.edu 3mm ^ 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 3mm

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.



§ 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


<cit.>

. In practice, the sample covariance[In 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


<cit.>

. See, e.g., <cit.> 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



<cit.>

. 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



<cit.>

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


<cit.>

. Spectral diagonal covariance models were successfully used in operational
statistical interpolation in meteorology in spherical geometry


<cit.>

, 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



<cit.>

.

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


<cit.>

, 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. <ref>, we provide a brief
statement of MLE and its asymptotic variance. In Sec. <ref>, 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. <ref>, and Sec. <ref>
contains computational illustrations. A comparison of the performance of MLE
for parametric models and of related shrinkage estimators is in
Sec. <ref>.







§ ASYMPTOTIC VARIANCE OF THE MAXIMUM LIKELIHOOD ESTIMATOR




First, we briefly review some standard results for reference. Suppose
𝕏_N=[  X_1,…,X_N]
is a random sample from a distribution on ℝ^n with density
f(  x,θ) with unknown parameter
vector θ in a parameter space Θ⊂ℝ
^p. The maximum likelihood estimate θ̂_N of the
true parameter θ^0 is defined by maximizing the
likelihood

    θ̂_N=max_θℒ(
    θ|𝕏_N)  ,  ℒ(
    θ|𝕏_N)  =∏_i=1^N
    ℒ(  θ|X_i)
    ,  ℒ(  θ|x)  =f(
    x,θ)  ,

or, equivalently, maximizing the log likelihood

    θ̂_N=max_θℓ(
    θ|𝕏_N)  ,  ℓ(
    θ|𝕏_N)  =∑_i=1^Nℓ(
    θ|X_i)  ,  ℓ(
    θ|x)  =log f(  x
    ,θ)  .



We adopt the usual assumptions that (i) the true parameter θ^0 lies in the interior of Θ, (ii) the density f determines the
parameter θ uniquely in the sense that f(x
,θ_1)=f(x,θ_2) a.s. if
and only if θ_1=θ_2, and (iii)
f(  x,θ) is a sufficiently
smooth function of x and θ. Then the error
of the estimate is asymptotically normal,

    √(N)(θ̂_N-θ^0
    )𝒩_p(0,Q_θ^0
    ), as N→∞,

where

    Q_θ^0=J_θ^0^-1,  
    J_θ^0=(  ∇_θ
    ℓ(θ^0|X)^⊤∇_θℓ(θ^0|X))  ,  X
    ∼ f(  x,θ^0)  .

The matrix J_θ^0 is called the Fisher information
matrix for the parameterization θ^0. Here,
X, x, and θ are columns,
while the gradient ∇_θℓ of ℓ with respect
to the parameter θ is a row vector, which is compatible
with the dimensioning of Jacobi matrices below. The mean value in
(<ref>) is taken with respect to X, which is
the only random quantity in (<ref>). Cf., e.g., <cit.> for details.



§ NESTED MAXIMUM LIKELIHOOD ESTIMATORS




Now suppose that we have an additional information that the true parameter
θ^0 lies in a subspace of Θ, which is
parameterized by k≤ p parameters (φ_1,…, φ
_k)^⊤=φ. Denote by ∇_φ
θ(φ) the p× k Jacobi matrix with
entries ∂θ_i/∂φ_j. In the next theorem,
we derive the asymptotic covariance of the maximum likelihood estimator for
φ,

    φ̂_N=max_φℓ(
    φ|𝕏_N)  ,  ℓ(
    φ|𝕏_N)  =∑_i=1^Nℓ(
    φ|X_i)  ,  ℓ(
    φ|x)  =log f(  x
    ,θ(  φ)  )  ,

based on the asymptotic covariance of θ in
(<ref>).


Assume that the map φ↦θ(φ) is one-to-one from Φ⊂ℝ^k to Θ, the map φ
↦θ(φ) is continuously
differentiable, ∇_φθ(φ) is full rank for all φ∈Φ,
and θ^0=θ(φ^0)
with φ^0 in the interior of Φ. Then,

    √(N)(φ̂_N-φ^0
    )𝒩_k(  0,Q_φ^0)   as N→∞,

where Q_φ^0=J_φ^0^-1, with
J_φ^0 the Fisher information matrix of the
parameterization φ given by

    J_φ^0=∇_φθ(φ^0)^⊤J_θ^0∇
    _φθ(φ^0).





From (<ref>) and the chain rule

    ∇_φℓ(φ|X
    )=∇_θℓ(θ|X
    )∇_φθ(φ),

we have

    J_φ^0     =(  ∇_φ
    ℓ(φ^0|X)^⊤∇
    _φℓ(φ^0|X)) 
         =∇_φθ(φ
    ^0)^⊤(  ∇_θℓ(θ
    ^0|X)^⊤∇_θℓ(θ^0|X))  ∇_φθ(φ^0)
         =∇_φθ(φ
    ^0)^⊤J_θ^0∇_φθ(φ^0).

The asymptotic distribution (<ref>) is now (<ref>)
applied to φ.


When the parameter θ is the quantity of interest in an
application, it is useful to express the estimate and its variance in terms of
the original parameter θ rather than the subspace
parameter φ.


Under the assumptions of Theorem <ref>,

    √(N)(θ(  φ̂_N)
    -θ^0)𝒩_p(  0
    ,Q_θ(  φ^0)  )
     as N→∞,

where

    Q_θ(  φ^0)  
    =∇_φθ(φ
    ^0)J_φ^0^-1∇_φ
    θ(φ^0)^⊤=∇
    _φθ(φ^0)(
    ∇_φθ(φ
    ^0)^⊤J_θ^0∇_φ
    θ(φ^0))  ^-1∇
    _φθ(φ^0)^⊤.





The lemma follows from (<ref>) by the delta method <cit.>, since the map φ↦θ
(φ) is continuously differentiable.



The matrix Q_θ(  φ^0)
is singular, so it cannot be written as the inverse of another matrix, but it
can be understood as the inverse J_θ(φ^0)^-1 of the Fisher information matrix for φ,
embedded in the larger parameter space Θ.


Suppose that ψ is another parameterization which satisfies
the same assumption as φ in Theorem <ref>:
the map ψ↦θ(ψ) is
one-to-one from Ψ⊂ℝ^m, k ≤ m ≤ p, to Θ,
ψ↦θ(ψ) is
continuously differentiable, ∇_ψθ(ψ) is full rank for all ψ∈Ψ, and
θ^0=θ(ψ^0), where
ψ^0 is in the interior of Ψ. Then, similarly as in
(<ref>), we have also

    √(N)(θ( ψ̂_N) -θ^0)𝒩_p(  0
    ,Q_θ(  ψ^0)  )  
    as N→∞,

where, as in (<ref>),

    Q_θ(  ψ^0)  =∇
    _ψθ(ψ^0
    )J_ψ^0^-1∇_ψθ(ψ^0)^⊤=∇_ψθ(ψ^0)(  ∇_ψθ(ψ^0)^⊤J_θ^0∇
    _ψθ(ψ^0))
    ^-1∇_ψθ(ψ
    ^0)^⊤.



The next theorem shows that when we have two parameterizations
φ and ψ which are nested, then the
smaller parameterization has smaller or equal asymptotic covariance than the
larger one. For symmetric matrices A and B, A≤ B means that A-B is
positive semidefinite.


Suppose that φ and
ψ satisfy the assumptions in Theorem <ref>,
and there exists a differentiable mapping φ↦ψ
from Φ to Ψ, such that φ^0↦ψ^0. Then,

    Q_θ(  φ^0)  ≤
    Q_θ(  ψ^0)  .

In addition, if U∼𝒩_p(  0
,Q_θ(  φ^0)  )
and V∼𝒩_p(  0,Q_θ(
ψ^0)  ) are random vectors with the
asymptotic distributions of the estimates θ(
φ̂_N) and θ(
ψ̂_N), then

    | U| ^2=1/N*Tr
    Q_θ(  φ^0)  ≤1/N*TrQ_θ(  ψ
    ^0)  =| V| ^2,

where | V| =(  V^⊤V)  ^1/2 is the
standard Euclidean norm in ℝ^p.



Denote A=J_θ^0, B=∇_φ
θ(φ^0), C=∇_ψ
θ(ψ^0). From the chain rule,

    ∇_φθ(  φ^0)  =∇_ψθ(
    ψ^0)  ∇_φψ(  φ^0)  ,

we have that B=C∇_φψ(
φ^0), and, consequently, *Range
B⊂*RangeC. Define

    P_B     =A^1/2B(B^⊤AB)^-1B^⊤A^1/2,
    
    P_C     =A^1/2C(C^⊤AC)^-1C^⊤A^1/2.

The matrices P_B and P_C are symmetric and idempotent, hence they are
orthogonal projections. In addition,

    *RangeP_B=*RangeA^1/2B⊂*RangeA^1/2C=*RangeP_C.

Consequently, P_B≤ P_C holds from standard properties of orthogonal
projections, and (<ref>) follows.

To prove (<ref>), note that for random vector X with
X=0 and finite second moment, | X| ^2
=*TrCovX from Karhunen-Loève decomposition
and Parseval identity. The proof is concluded by using the fact that for
symmetric matrices, A≤ B implies *TrA≤*TrB, cf. e.g., <cit.>.



In the practically interesting cases when there is a large difference in the
dimensions of the parameters φ and ψ,
many eigenvalues in the covariance of the estimation error become zero. The
computational tests in Sec. <ref> show that the resulting decrease of
the estimation error can be significant.




§ 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 <cit.> and
Gaussian Markov Random Fields (GMRF), <cit.>, 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


<cit.>

. The results of Section <ref> apply for this case and we shall
illustrate them in the simulation study of Section <ref>.

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
<cit.>. 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 <cit.> and we shall treat
the diagonal case in more detail.

Suppose

    X∼𝒩_n(0,D),

where 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


<cit.>

, 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.



 §.§ Sample covariance




The top-level parameter space Θ consists of all symmetric positive
definite matrices, resulting in the parameterization Σ with
n(  n+1)  /2 independent parameters. The likelihood of a
sample 𝕏_N=[  X^(1),…,X
^(N)] from 𝒩_n(0,Σ) is

    L(  Σ|𝕏_N)  =1/(  Σ)
    ^N/2(  2π)  ^nN/2e^-1/2*Tr(
    Σ^-1𝕏_N𝕏_N^⊤)  .

If N≥ n, it is well known (e.g. <cit.>, p. 83) that the
likelihood is maximized at what we call here sample covariance matrix

    Σ̂_N=1/N∑_i=1^NX^(i)(
    X^(i))  ^⊤.

The Fisher information matrix of the sample covariance estimator is
<cit.>

    J^(0)((Σ))=1/2Σ^-1⊗Σ^-1,

where ⊗ stands for the Kronecker product and  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× n^2.


 If Σ̂_N is singular, L(  Σ̂_N|𝕏_N) cannot be evaluated because that requires the
inverse of Σ̂. Also, in this case the likelihood L(
Σ|𝕏_N) is not bounded above on the set of all
Σ>0, thus the maximum of L(  Σ|𝕏_N) does
not exist on that space. To show that, consider an orthonormal change of basis
so that the vectors in *span(  𝕏_N)
come first, write vectors and matrices in the corresponding 2×2 block
form, and let

    Σ̃_N=[
    [ Σ̃_11     0;     0     0 ]]  ,  Σ̃_11>0.

Then lim_a→0^+𝕏_N^⊤(  Σ̃
_N+aI)  ^-1𝕏_N exists, but lim_a→0^+
(  Σ̃_N+aI)  =0, thus

    lim_a→0^+L(  Σ̃_N+aI|𝕏_N)=∞.



Note that when the likelihood is redefined in terms of the subspace
*span(  𝕏_N) only, the sample
covariance can be obtained by maximization on the subspace <cit.>.


When the true covariance is diagonal (Σ≡ D, cf. (<ref>)), a
significant improvement can be achieved by setting the off-diagonal terms of
sample covariance to zero,

    D̂_N^(  0)  =diag(  Σ̂
    _N)  .

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

    Σ̂_N^u=1/N-1∑_i=1^NX^(i)(
    X^(i))  ^⊤

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

    |D̂_N^(  0)  -D| _F≤|Σ̂_N^u-D| _F

cf. <cit.> for the case when the mean is assumed to be known
like here, and <cit.> for the unbiased sample covariance
and unknown mean.



 §.§ Diagonal covariance




The parameter space Θ_1 consisting of all diagonal matrices with
positive diagonal, with n parameters d=(  d_1
,…,d_n)  ^⊤, can be viewed as a simple class of models for
either covariance or its inverse. The log-likelihood function for
D=(d_1,…,d_n) with a given random sample 𝕏
_N=[  X^(1),…,X^(N)] from
𝒩_n(  0,D) is

    ℓ(D|𝕏_N)=-N/2log(  (2π)^n|D|)  -1/2∑_k=1^N(  X^(k))  ^⊤
    D^-1X^(k)

and has its maximum at

    d̂_j=1/N∑_k=1^N(  X_j^(k))  ^2
    ,j=1,…,n,

where X_j^(k) denotes the j-th entry of X^(k). The
sum of squares S_j^2=∑_k=1^N(  X_j^(k))  ^2 is a
sufficient statistic for the variance d_j. Thus, we get the maximum
likelihood estimator

    D̂_N^(  1)  =1/N(  S_1^2
    ,…,S_n^2)  .



It is easy to compute the Fisher information matrix explicitly,

    J_D^(  1)  =(  1/2d_1^2,…,1/2d_n^2)  .

which is an n × n matrix and gives the asymptotic covariance of the
estimation error

    1/NQ_D^(  1)  =1/NJ_D^(  1)  
    ^-1=1/N(  2d_1^2,…, 2d_n^2)

from (<ref>).



 §.§ 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_2h_i)f_i(α))^-1, i=1,…,n, where
c_1,c_2 and α are unknown parameters, h_1,…,h_n are
known positive numbers, and f_1,…,f_n are known differentiable
functions. For easier computation it is useful to work with τ_i=1/d_i=(c_1+c_2h_i)f_i(α). Maximum likelihood estimators for
c_1,c_2, and α can be computed effectively from the likelihood

    ℓ(D|𝕏_N)=-N/2nlog(2π)+N/2∑_i=1^n
    logτ_i-1/2∑_i=1^nτ_iS_i^2

by using the chain rule. It holds that

    ∂ℓ/∂ c_1     =∑_i=1^n∂ℓ/∂τ_i∂τ_i/∂ c_1=∑_i=1
    ^n(  N/2τ_i-S_i^2/2)  ∂τ_i/∂ c_1
         =N/2∑_i=1^n(  1/(c_1+c_2h_i)f_i(α
    )-1/NS_i^2)  f_i(α).

Setting this derivative equal to zero we get

    ∑_i=1^n(  1/c_1+c_2h_i-1/NS_i^2
    f_i(α))  =0.

Analogously,

    ∂ℓ/∂ c_2=∑_i=1^n∂ℓ/∂τ_i∂τ_i/∂ c_2=N/2
    ∑_i=1^n(  1/(c_1+c_2h_i)f_i(α)-1/NS_i^2)  h_if_i(α),

so the equation for estimating the parameter c_2 is

    ∑_i=1^n(  h_i/c_1+c_2h_i-1/NS_i^2
    h_if_i(α))  =0.

Similarly,

    ∂ℓ/∂α     =∑_i=1^n∂ℓ/∂τ_i∂τ_i/∂α=N/2
    ∑_i=1^n(  1/(c_1+c_2h_i)f_i(α)-1/NS_i^2)  (c_1+c_2h_i)∂ f_i(α)/∂α
         =N/2∑_i=1^n(  1/f_i(α)-1/N
    S_i^2(c_1+c_2h_i))  ∂ f_i(α)/∂α

and setting the derivative to zero, we get

    ∑_i=1^n(  1/f_i(α)∂ f_i(α
    )/∂α-1/NS_i^2(c_1+c_2h_i)∂
    f_i(α)/∂α)  =0.



The maximum likelihood estimator for D is then given by

    D̂^(3)={((ĉ_1+ĉ_2h_i)f_i(α̂))^-1,
    i = 1,…,n },

where (ĉ_1,ĉ_2,α̂) is the solution of the system
(<ref>), (<ref>), (<ref>). This expression corresponds to
searching a maximum likelihood estimator of D in the subspace Θ
_3⊂Θ_1⊂Θ formed by diagonal matrices
{((c_1+c_2h_i)f_i(α))^-1,i=1,…,n}.

For completeness, the asymptotic covariance of the estimation error about

    D^(3)={d_i(c_1,c_2,α), i=1,…,n},

contained in 𝕏_N is

    1/NQ_D^(  3)  =1/N∇d
    (c_1,c_2,α)J_c_1,c_2,α^-1∇d
    (c_1,c_2,α)^⊤

from (<ref>), where the Fisher information matrix J_c_1
,c_2,α is the 3×3 matrix

    J_c_1,c_2,α=
       [                    1/2∑_i=1^n1/(c_1+c_2h_i)^2                 1/2
    ∑_i=1^nh_i/(c_1+c_2h_i)^2   1/2∑_i=1
    ^n1/(c_1+c_2h_i)f_i(α)∂ f_i(α
    )/∂α;                  1/2∑_i=1^nh_i/(c_1+c_2h_i)^2                1/2∑_i=1^nh_i^2/(c_1+c_2h_i)^2  1/2
    ∑_i=1^nh_i/(c_1+c_2h_i)f_i(α)∂
    f_i(α)/∂α;    1/2∑_i=1^n1/(c_1+c_2h_i)f_i(α)
    ∂ f_i(α)/∂α 1/2∑_i=1^n
    h_i/(c_1+c_2h_i)f_i(α)∂ f_i(α
    )/∂α      1/2∑_i=1^n1/f_i^2(α
    )(  ∂ f_i(α)/∂α)  ^2 ]

and

    d(c_1,c_2,α)       =[  d_1(c_1,c_2
    ,α),…,d_n(c_1,c_2,α)]  ^⊤
         =[  ((c_1+c_2h_1)f_1(α))^-1,…,((c_1+c_2
    h_n)f_n(α))^-1]  ^⊤.





 §.§ 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(α))^-1, i.e. τ_i=cf_i(α),
i=1,…,n, where c and α are unknown parameters. Maximum
likelihood estimators for c and α can be computed similarly as in the
previous case. The estimating equations have the form

    1/c     =1/n∑_i=1^n1/NS_i^2f_i(α)
    1/c∑_i=1^n1/f_i(α)∂ f_i(α
    )/∂α     =∑_i=1^n1/NS_i^2∂
    f_i(α)/∂α,

which can be rearranged to

    1/c     =1/n∑_i=1^n1/NS_i^2f_i
    (α)
    
    0       =∑_i=1^nS_i^2f_i(α)(  1/f_i(α
    )∂ f_i(α)/∂α-1/n∑_j=1^n
    1/f_j(α)∂ f_j(α)/∂α)  .

Equation (<ref>) is an implicit formula for estimating α. Its
result can be used for estimating c through (<ref>). The maximum
likelihood estimator for D is then given by

    D̂^(2)=(  (ĉf_1(α̂))^-1,…,(ĉf_n(α̂))^-1)  ,

where ĉ and α̂ are MLEs of c and α. It
corresponds to searching a maximum likelihood estimator of D in the subspace
Θ_2⊂Θ_3⊂Θ_1⊂Θ formed by diagonal
matrices {  (cf_i(α))^-1,i=1,…,n }. Of
course, the estimator D̂^(2) does not have “larger" variance than
D̂^(3).

The covariance of the asymptotic distribution of the parameters d_1
,…,d_n is

    1/NQ_D^(  2)  =1/N∇d
    (c,α)J_c,α^-1∇d(c,α)^⊤,

from (<ref>), where Fisher information matrix at D={d_i
(c,α),i=1,…,n} is the 2×2 matrix

    J_c,α=
    [                                   n/2c^2          1/2c∑_i=1^n1/f_i(α)
    ∂ f_i(α)/∂α;          1/2c∑_i=1^n1/f_i(α)∂ f_i(α
    )/∂α 1/2∑_i=1^n1/f_i^2(α
    )(  ∂ f_i(α)/∂α)  ^2 ]

and d(c,α)=[  d_1(c,α),…,d_n
(c,α)]  ^⊤=[  (cf_1(α))^-1,…,(cf_n
(α))^-1]  ^⊤.



 §.§ 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
<cit.> and applied on data from oceanography. The corresponding
Fisher information matrix may be found as the negative of the Hessian matrix
<cit.>.






§ COMPUTATIONAL STUDY




In Section <ref>, 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.






 §.§ 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×10 (resulting in n=100) and inverse covariance structure as in
Fig. <ref>. 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 Σ^-1 with 4 neighbors was set as
“truth” and random samples were generated
from 𝒩_n(0,Σ) with sample sizes
N=10,15,20,…,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 <cit.>.

The difference of each estimator from the true matrix Σ 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. <ref>.






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.



 §.§ 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=1/ce^αλ_i,i=1,…,n,
where c and α are parameters and λ_i are the eigenvalues of
Laplace operator in two dimensions on 10×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 𝒩_n(0,D)
with sample sizes N=5,…,20. For each sample, four covariance matrix
estimators were computed:



  * sample covariance matrix Σ̂_N, cf.
(<ref>)


  * diagonal part D̂^(  0) of the sample covariance
matrix, cf. (<ref>)


  * MLE D̂^(  1) in the space of diagonal matrices,
cf. (<ref>)


  * MLE D̂^(  3)  ={(ĉ_1-ĉ_2
λ_i)^-1e^α̂λ_i,i=1,…,n} with 3 parameters
c_1,c_2 and α, cf. (<ref>).


  * MLE D̂^(  2)  ={ĉ^-1e^α̂λ_i,i=1,…,n} with 2 parameters c and α, cf.
(<ref>).


Let us briefly discuss the choice of the covariance model d_i=1/ce^αλ_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.,
<cit.>). Exponential decay is used, e.g., in
<cit.>. 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 -λ_i,i=1,…,n, e.g., <cit.>.

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. <ref>.









For the diagonal MLE, given by (<ref>), (<ref>), and (<ref>), we
can expect from (<ref>) that these estimators should
satisfy asymptotically

    (  |D̂_N^(  k)  -D| _F
    ^2)  ≈1/N*Tr(J_D^(  k)  
    ^-1),   k=1,2,3,

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

    *Tr(J_D^(  2)  ^-1)≤*Tr
    (J_D^(  3)  ^-1)≤*Tr(J_D^(
    1)  ^-1)

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

    |D̂_N^(  2)  -D| _F^2
    ≤|D̂_N^(  3)  -D| _F^2
    ≤|D̂_N^(  1)  -D| _F^2,

which is confirmed by the simulations (see Figure <ref>, resp.
<ref>).

The comparisons (<ref>) of the Frobenius norm of the
error in the mean squared complement the pointwise comparison
(<ref>) 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, Σ̂_N is singular and
cannot be cast as MLE with an accompanying Fisher information matrix, cf.
Remark <ref>. 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.









§ 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 Σ̂_N and the so called target matrix T:

    Σ̂_S=γΣ̂_N+(1-γ)T, for
    γ∈0,1].





One of the simplest shrinkage estimators has the form of (<ref>
) with the target matrix equal to identity, which results in shrinking all
sample eigenvalues with the same intensity towards their mean value.
<cit.> derived the optimal shrinkage parameter γ to
minimize the squared Frobenius loss

    min_γ||Σ̂_S-D||_F^2.

The comparison of this estimator with the maximum likelihood estimator
D̂^(2) was accomplished by a simulation with identical setting as in
Section <ref>. The results are shown in Fig. <ref>.
For reference, the sample covariance Σ̂_N and its diagonal
D̂^(0) are also added.

Another regularization method is described in

<cit.>. 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 κ_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 κ_max, i.e.

    max_Σℓ(Σ) subject to λ
    _max(Σ)/λ_min(Σ)≤κ_max,

where λ_max(Σ), resp. λ_min(Σ), is the largest,
resp. the smallest, eigenvalue of the covariance matrix Σ. An optimal
κ_max is selected by maximization of the expected likelihood, which
is approximated by using K-fold cross-validation. The authors proved that
κ_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 <cit.>. In
Fig. <ref>, we can see the performance of this estimator
(denoted as ) in comparison of other methods.




The shrinkage estimator Σ̂_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
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 Σ̂_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. <ref> that the
condition-number-regularized estimator provides more precise estimates than
the sample covariance Σ̂_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
κ_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.
















































































§ 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.


apacite


