эффективный способ инвертировать большую корреляционную матрицу

Допустим, у меня есть очень большая корреляционная матрица такой формы:

t1.rep1 = rnorm(n=100,mean=10,sd=)
t2.rep1 = t1.rep1 + rnorm(n=100,mean=3,sd=2)
t3.rep1 = t1.rep1 + rnorm(n=100,mean=2,sd=2)
t1.rep2 = rnorm(n=100,mean=2,sd=1)
t2.rep2 = t1.rep2 + rnorm(n=100,mean=0.5,sd=0.5)
t3.rep2 = t1.rep2 + rnorm(n=100,mean=0.7,sd=0.9)
t1.rep3 = rnorm(n=100,mean=2,sd=1)
t2.rep3 = t1.rep3 + rnorm(n=100,mean=0.5,sd=0.5)
t3.rep3 = t1.rep3 + rnorm(n=100,mean=0.7,sd=0.9)
Sigma = matrix(
  c(cov(t1.rep1, t1.rep1), 0, 0, cov(t1.rep1, t2.rep1), 0, 0, cov(t1.rep1, t3.rep1), 0, 0,
  0, cov(t1.rep2, t1.rep2), 0, 0, cov(t1.rep2, t2.rep2), 0, 0, cov(t1.rep2, t3.rep2), 0,
  0, 0, cov(t1.rep3, t1.rep3), 0, 0, cov(t1.rep3, t2.rep3), 0, 0, cov(t1.rep3, t3.rep3),
  cov(t2.rep1, t1.rep1), 0, 0, cov(t2.rep1, t2.rep1), 0, 0, cov(t2.rep1, t3.rep1), 0, 0,
  0, cov(t2.rep2, t1.rep2), 0, 0, cov(t2.rep2, t2.rep2), 0, 0, cov(t2.rep2, t3.rep2), 0,
  0, 0, cov(t2.rep3, t1.rep3), 0, 0, cov(t2.rep3, t2.rep3), 0, 0, cov(t2.rep3, t3.rep3),
  cov(t3.rep1, t1.rep1), 0, 0, cov(t3.rep1, t2.rep1), 0, 0, cov(t3.rep1, t3.rep1), 0, 0,
  0, cov(t3.rep2, t1.rep2), 0, 0, cov(t3.rep2, t2.rep2), 0, 0, cov(t3.rep2, t3.rep2), 0,
  0, 0, cov(t3.rep3, t1.rep3), 0, 0, cov(t3.rep3, t2.rep3), 0, 0, cov(t3.rep3, t3.rep3)),
  nrow = 9, ncol = 9)

Моя корреляционная матрица — сигма.

И я хочу вычислить его инверсию, т.е.

Sigma.inv = solve(Sigma)

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

Есть ли способ использовать разреженность и структуру матрицы для вычисления обратной сигмы гораздо более быстрым/эффективным способом?

Мне нужно вычислить это обратное итеративно, поэтому быстрый способ вычисления обратного значительно поможет моему алгоритму ускориться.


r covariance inverse
person Dnaiel    schedule 25.02.2014    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Предоставленный вами Sigma на самом деле является блочно-диагональным:

x = c(1,4,7,2,5,8,3,6,9)
Sigma[x,x]

          [,1]     [,2]      [,3]     [,4]     [,5]     [,6] ...
[1,] 0.9494388 1.130673 0.9825316 0.000000 0.000000 0.000000 ...
[2,] 1.1306727 4.983144 1.2112634 0.000000 0.000000 0.000000 ...
[3,] 0.9825316 1.211263 5.0771423 0.000000 0.000000 0.000000 ...
[4,] 0.0000000 0.000000 0.0000000 1.211892 1.223293 1.328587 ...
[5,] 0.0000000 0.000000 0.0000000 1.223293 1.469146 1.242400 ...
[6,] 0.0000000 0.000000 0.0000000 1.328587 1.242400 2.377406 ...
...

А блочно-диагональные матрицы можно инвертировать намного быстрее поблочно. Просто замените каждый блок его инверсией.

person Andrey Shabalin    schedule 25.02.2014
comment
@Шабалин, большое спасибо, похоже, отличное предложение. Есть ли какой-нибудь трюк, чтобы быстро сделать инверсию каждого блока в R? Просто применить функцию применения к каждому блоку? или сделать цикл for? for циклы, как правило, довольно медленные в R... - person Dnaiel; 25.02.2014
comment
Может быть, вы можете хранить матрицу в блоках для начала. - person Andrey Shabalin; 26.02.2014
comment
Кроме того, взгляните на stats.stackexchange.com/questions/14951/ - person Andrey Shabalin; 26.02.2014
comment
@Shabalin Большое спасибо, все указатели в правильном направлении. Однако поблочная инверсия все еще довольно медленная. Я попробовал холески, и это было не очень хорошо. Есть ли способ еще больше ускорить инверсию матрицы, хорошо ли работает многопоточность? Я опубликую это как отдельный вопрос. - person Dnaiel; 26.02.2014