Невозможно получить положительно определенную матрицу дисперсии при очень малых собственных значениях

Для запуска канонического анализа соответствия (cca package ade4) мне нужна положительно определенная дисперсионная матрица. (Что в теории всегда так), но:

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

Последнее собственное значение равно -1,377693e-09, что равно ‹ 0. Но теоретическое значение > 0.
Я не могу запустить функцию, если одно из собственных значений равно ‹ 0.

Я действительно не знаю, как это исправить, не меняя код функции cca()

Спасибо за помощь


r matrix statistics eigenvalue
person joel lapalme    schedule 31.05.2013    source источник
comment
***На самом деле теоретическое значение всегда ›= 0 ***   -  person joel lapalme    schedule 31.05.2013

Ответы (2)


arrow_upward
3
arrow_downward

Вот два подхода:

V <- var(a)

# 1
library(Matrix)
nearPD(V)$mat

# 2 perturb diagonals
eps <- 0.01
V + eps * diag(ncol(V))
person G. Grothendieck    schedule 31.05.2013

arrow_upward
4
arrow_downward

Вы можете немного изменить ввод, чтобы сделать матрицу положительно определенной.

Если у вас есть матрица дисперсии, вы можете обрезать собственные значения:

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

Используя разложение по сингулярным числам, вы можете сделать то же самое непосредственно с 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

Однако это изменяет вашу матрицу a на 0.1, что очень много (я подозреваю, что это так много, потому что матрица a квадратная: в результате одно из собственных значений var(a) равно 0).

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

На самом деле мы можем добиться большего, просто добавив немного шума.

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
person Vincent Zoonekynd    schedule 31.05.2013