# How to avoid 0 determinant when sample covariance matrix has very small values

I have $n$, $p$-dimensional vectors and I am constructing the $p \times p$ covariance matrix using the following formula:

$$\mathrm{Cov}(j,k) = \frac{1}{n-1} {\sum^n_{i=1}} (x_i(j) - {\mu}(j))(x_i(k) - {\mu(k)})$$

However, the values in the $p$-dimensional vectors are between 0 and 1, as it is a normalized vector, and $n$ can be large which results in $\mathrm{Cov}(j,k)$ having rather small elements.

I am trying to calculate the determinant of the covariance matrix to determine the differential entropy of the underlying multivariate Gaussian. However, because of the small numbers, the determinant is effectively 0 (and in fact rounds to 0 in computation).

Are there any methods to avoid this? Perhaps I can avoid the division by $(n-1)$?

Several things.

• The determinant is scale equivariant, so without loss of generality you can rescale your $x_i$'s (say multiply them by 10) and then factor back the scale after all the intermediate operations have been performed as shown in the following R code:

p <- 10
n <- 100
a1 <- matrix(runif(n*p),n,p)
a1 <- sweep(a1,1,sqrt(rowSums(a1*a1)),FUN="/")
a2 <- a1*10 #you can re-scale your data
log(det(cov(a1[1:(p+1),])))
log(det(cov(a2[1:(p+1),])))-log(10^(2*(p))) #and compensate for the rescaling 'at the end'

• You should work with the determinant rescaled to be expressed in the original units of your dataset (e.g. if $\varSigma$ is your covariance matrix and $p$ it's rank, use $|\varSigma|^{1/2p}$). This is similar to how the standard deviations are expressed in the original units but the variances in the square of the original units.

• What is matrix(runif(n*p),n,p)?
– Aly
Feb 19, 2013 at 11:28
• @Aly: your data (a bunch of random points between 0 and 1). The next line normalizes it... Feb 19, 2013 at 11:40
• so essentially $log(det(cov)) == log((1/{10^2})*det(cov*10)) == log(det(cov*10)) - 2log(10)$ ?
– Aly
Feb 19, 2013 at 11:45
• @Aly: you forgot the p! Feb 19, 2013 at 11:46
• Ah yes, I almost forgot what I was trying to do and just focusing on the equation. I understand now
– Aly
Feb 19, 2013 at 12:09

Numerical stability will be a problem for you.

One approach to make the calculation a little more stable - taking $X$ to be a mean-corrected matrix of the $x$-values, I think it makes sense to either form a QR decomposition of the $X$-matrix or failing that a Choleski decomposition of the $X'X$ matrix, and use the product of the diagonal of $R$ (or of the Choleski factor) which will be the square root of the determinant.

There are other approaches (some of which may be better than computing the QR).

You still need to scale for the $n-1$'s, though.