It is not possible to obtain a positive definite dispersion matrix when very small eigenvalues

Question

To perform canonical conformance analysis (cca package ade4), I need a positive definite variance matrix. (Which in theory always holds) but:

matrix(c(2,59,4,7,10,0,7,0,0,0,475,18714,4070,97,298,0,1,0,17,7,4,1,4,18,36),nrow=5) > a [,1] [,2] [,3] [,4] [,5] [1,] 2 0 475 0 4 [2,] 59 7 18714 1 1 [3,] 4 0 4070 0 4 [4,] 7 0 97 17 18 [5,] 10 0 298 7 36 > eigen(var(a)) $values [1] 6.380066e+07 1.973658e+02 3.551492e+01 1.033096e+01 [5] -1.377693e-09 

The last eigenvalue is -1.377693e-09 , which is <0. But the theoretical value is → 0.
I cannot run the function if one of the eigenvalues ​​is <0

I really don't know how to fix this without changing the cca () function code

thanks for the help

Answer 1

Here are two approaches:

 V <- var(a) # 1 library(Matrix) nearPD(V)$mat # 2 perturb diagonals eps <- 0.01 V + eps * diag(ncol(V)) 

Answer 2

You can change the input a little to make the matrix positive.

If you have a variance matrix, you can trim your eigenvalues:

 correct_variance <- function(V, minimum_eigenvalue = 0) { V <- ( V + t(V) ) / 2 e <- eigen(V) e$vectors %*% diag(pmax(minimum_eigenvalue,e$values)) %*% t(e$vectors) } v <- correct_variance( var(a) ) eigen(v)$values # [1] 6.380066e+07 1.973658e+02 3.551492e+01 1.033096e+01 1.326768e-08 

Using the expansion of singular values, you can do the same directly with a .

 truncate_singular_values <- function(a, minimum = 0) { s <- svd(a) s$u %*% diag( ifelse( s$d > minimum, s$d, minimum ) ) %*% t(s$v) } svd(a)$d # [1] 1.916001e+04 4.435562e+01 1.196984e+01 8.822299e+00 1.035624e-01 eigen(var( truncate_singular_values(a,.2) ))$values # [1] 6.380066e+07 1.973680e+02 3.551494e+01 1.033452e+01 6.079487e-09 

However, this changes your matrix a by 0.1 , which is a lot (I suspect it is high because the matrix a is square: as a result, one of the eigenvalues ​​of var(a) exactly 0.)

 b <- truncate_singular_values(a,.2) max( abs(ba) ) # [1] 0.09410187 

In fact, we can do better by just adding some noise.

 b <- a + 1e-6*runif(length(a),-1,1) # Repeat if needed eigen(var(b))$values # [1] 6.380066e+07 1.973658e+02 3.551492e+01 1.033096e+01 2.492604e-09