компактная / эффективная замена для diag (XVX ^ T)? - PullRequest
Новая
       25

компактная / эффективная замена для diag (XVX ^ T)?

6 голосов
/ 10 июля 2019

При прогнозировании линейной статистической модели у нас обычно есть матрица модели X предикторов, соответствующих точкам, в которых мы хотим делать прогнозы; вектор коэффициентов beta; и матрица дисперсии-ковариации V. Вычисление прогнозов просто X %*% beta. Самый простой способ вычисления дисперсий прогнозов - это 

diag(X %*% V %*% t(X))

или чуть эффективнее 

diag(X %*% tcrossprod(V,X))

Однако, это очень неэффективно, потому что оно строит матрицу n * n, когда все, что мы действительно хотим, это диагональ. Я знаю, что мог бы написать какую-нибудь Rcpp-петлю, которая бы вычисляла только диагональные члены, но мне интересно, есть ли существующий трюк с линейной алгеброй в R, который будет хорошо делать то, что я хочу ... (если кто-то хочет написать Rcpp-loopy для меня в качестве ответа я бы не стал возражать, но я бы предпочел решение на основе чистого R)

FWIW predict.lm, кажется, делает что-то умное, умножая X на инверсию R-компонента QR-разложения lm; Я не уверен, что это всегда будет доступно, но это может быть хорошей отправной точкой (см. здесь )

Ответы [ 2 ]

3 голосов
/ 10 июля 2019

В соответствии с этим вопросом Octave / Matlab , для двух матриц A и B мы можем использовать тот факт, что диагональная запись nth для AB будетпроизведение строки nth на A со столбцом nth на B.Мы можем наивно распространить это на случай трех матриц ABC.Я не рассматривал, как оптимизировать в случае, когда C=A^T, но, кроме этого, этот код выглядит как многообещающее ускорение: 

start_time <- Sys.time()

A=matrix(1:1000000, nrow = 1000, ncol = 1000)
B=matrix(1000000:1, nrow = 1000, ncol = 1000)

# Try one of these two
res=diag(A %*% B %*% t(A)) # ~0.47s
res=rowSums(A * t(B %*% t(A))) # ~0.27s

end_time <- Sys.time()

print(end_time - start_time)

Использование tcrossprod не ускорило результаты при запускеэтот кодТем не менее, просто использование подхода row-sum-dot-product кажется гораздо более эффективным, по крайней мере, на этом глупом примере, который предполагает (хотя я не уверен), что rowSums не вычисление полных промежуточных матриц перед возвратом диагональных элементов, как я и ожидал, происходит с diag.

0 голосов
/ 10 июля 2019

Я не совсем уверен, насколько это эффективно,

  1. Найдите U такой, что V = U %*% t(U); это возможно, поскольку V - матрица ков.
  2. XU = X %*% U
  3. result = apply(XU, 1, function(x) sum(x^2))

Демо 

V <- cov(iris[, -5])
X <- as.matrix(iris[1:5, -5])

Использование SVD 

svd_v <- svd(V)
U <- svd_v$u %*% diag(sqrt(svd_v$d))
XU = X %*% U
apply(XU, 1, function(x) sum(x^2))
#       1        2        3        4        5 
#41.35342 39.36286 35.42369 38.25584 40.30839

Другой подход - это также не будет быстрее, чем у @ davewy's 

U <- chol(V)
XU = (X %*% U)^2
rowSums(XU)
...