How to ensure properties of covariance matrix when fitting multivariate normal model using maximum likelihood?

Question

Suppose I have the following model

$$y_i=f(x_i,\theta)+\varepsilon_i$$

where $y_i\in \mathbb{R}^K$ , $x_i$ is a vector of explanatory variables, $\theta$ is the parameters of non-linear function $f$ and $\varepsilon_i\sim N(0,\Sigma)$, where $\Sigma$ naturally is $K\times K$ matrix.

The goal is the usual to estimate $\theta$ and $\Sigma$. The obvious choice is maximum likelihood method. Log-likelihood for this model (assuming we have a sample $(y_i,x_i),i=1,...,n$) looks like

$$l(\theta,\Sigma)=-\frac{n}{2}\log(2\pi)-\frac{n}{2} \log\det\Sigma-\sum_{i=1}^n(y_i-f(x_i,\theta))'\Sigma^{-1}(y-f(x_i,\theta)))$$

Now this seems simple, the log-likelihood is specified, put in data, and use some algorithm for non-linear optimisation. The problem is how to ensure that $\Sigma$ is positive definite. Using for example optim in R (or any other non-linear optimisation algorithm) will not guarantee me that $\Sigma$ is positive definite.

So the question is how to ensure that $\Sigma$ stays positive definite? I see two possible solutions:

1. Reparametrise $\Sigma$ as $RR'$ where $R$ is upper-triangular or symmetric matrix. Then $\Sigma$ will always be positive-definite and $R$ can be unconstrained.

2. Use profile likelihood. Derive the formulas for $\hat\theta(\Sigma)$ and $\hat{\Sigma}(\theta)$. Start with some $\theta_0$ and iterate $\hat{\Sigma}_j=\hat\Sigma(\hat\theta_{j-1})$, $\hat{\theta}_j=\hat\theta(\hat\Sigma_{j-1})$ until convergence.

Is there some other way and what about these 2 approaches, will they work, are they standard? This seems pretty standard problem, but quick search did not give me any pointers. I know that Bayesian estimation would be also possible, but for the moment I would not want to engage in it.

Answer 1

As it turns out you can use profile maximum likelihood to ensure the necessary properties. You can prove that for given $\hat\theta$, $l(\hat\theta,\Sigma)$ is maximised by

$$\hat\Sigma=\frac{1}{n}\sum_{i=1}^n\hat{\varepsilon}_i\hat{\varepsilon}_i',$$

where

$$\hat{\varepsilon}_i=y_i-f(x_i,\hat\theta)$$

Then it is possible to show that

$$\sum_{i=1}^n(y_i-f(x_i,\hat\theta))'\hat\Sigma^{-1}(y-f(x_i,\hat\theta)))=const,$$

hence we only need to maximise

$$l_R(\theta,\Sigma)=-\frac{n}{2} \log\det\hat\Sigma.$$

Naturally in this case $\Sigma$ will satisfy all the necessary properties. The proofs are identical for the case when $f$ is linear which can be found in Time Series Analysis by J. D. Hamilton page 295, hence I omitted them.

Answer 2

Assuming that in constructing the covariance matrix, you are automatically taking care of the symmetry issue, your log-likelihood will be $-\infty$ when $\Sigma$ is not positive definite because of the $\log {\rm det} \ \Sigma$ term in the model right? To prevent a numerical error if ${\rm det} \ \Sigma < 0$ I would precalculate ${\rm det} \ \Sigma$ and, if it is not positive, then make the log likelihood equal -Inf, otherwise continue. You have to calculate the determinant anyways, so this is not costing you any extra calculation.

Answer 3

An alternative parameterization for the covariance matrix is in terms of eigenvalues $\lambda_1,...,\lambda_p$ and $p(p-1)/2$ "Givens" angles $\theta_ij$.

That is, we can write

$$\Sigma = G^T \Lambda G$$

where $G$ is orthonormal, and

$$\Lambda = diag(\lambda_1, ..., \lambda_p)$$

with $\lambda_1 \geq ... \geq \lambda_p \geq 0$.

Meanwhile, $G$ can be parameterized uniquely in terms of $p(p-1)/2$ angles, $\theta_{ij}$, where $i = 1,2,...,p-1$ and $j = i, ..., p-1$.[1]

(details to be added)

[1]: Hoffman, Raffenetti, Ruedenberg. "Generalization of Euler Angles to N‐Dimensional Orthogonal Matrices". J. Math. Phys. 13, 528 (1972)

Answer 4

Along the lines of charles.y.zheng's solution, you may wish to model $\Sigma = \Lambda + C C^{\top}$, where $\Lambda$ is a diagonal matrix, and $C$ is a Cholesky factorization of a rank update to $\Lambda$. You only then need to keep the diagonal of $\Lambda$ positive to keep $\Sigma$ positive definite. That is, you should estimate the diagonal of $\Lambda$ and the elements of $C$ instead of estimating $\Sigma$.