Блог

Главная - Вопросы по программированию - R: Запись двойного суммирования в матричном виде

R: Запись двойного суммирования в матричном виде

#r #matrix #vectorization #linear-algebra

Вопрос:

введите описание изображения здесь

Я хотел бы избавиться от суммирования по $t$ и реализовать это в R. До сих пор у меня есть это

 Z_diff = mapply(FUN = function(i) {Z-Z[i]}, i=1:n, SIMPLIFY = TRUE)
rss = 0
for(t in 1:n){
    rss = rss   t(Y - X %*% B[,t]) %*% (diag(c(Z_diff[,t])) %*% (Y - X %*% B[,t]))
  }

Комментарии:

1. И ваш вопрос в чем?

2. @deschen Как я могу избавиться от суммирования более $t$, чтобы избежать использования for цикла?

Ответ №1:

Ниже мы определяем тестовый ввод и вычисляем rss, используя исходную формулу в вопросе. Затем мы показываем матричную формулу, которая дает тот же результат. Если вы выпишете формулы для компонентов матрицы в сумме последней строки, причина результата должна быть ясна.

 # input
set.seed(123)
n <- 4
Z <- rnorm(n)
Y <- rnorm(n)
B <- matrix(rnorm(n*n), n)
X <- matrix(rnorm(n*n), n)

# scalar sum
rss <- 0
for(i in 1:n) 
  for(t in 1:n) 
    rss <- rss   (Y[i] - c( crossprod(X[i, ], B[, t])))^2 * (Z[t] - Z[i])
rss
## [1] 43.39269

# matrix sum
ZZ <- t(outer(Z, Z, "-"))
sum(ZZ * (Y - X%*%B)^2)
## [1] 43.39269

Производительность

С приведенными выше данными формулировка матрицы на моей машине работает примерно в 300 раз быстрее.

 library(microbenchmark)

microbenchmark(
scalar = {
  rss <- 0
  for(i in 1:n) 
    for(t in 1:n) 
      rss <- rss   (Y[i] - c( crossprod(X[i, ], B[, t])))^2 * (Z[t] - Z[i])
  rss
},
matrix = {
  ZZ <- t(outer(Z, Z, "-"))
  sum(ZZ * (Y - X%*%B)^2)
})
  })
## Unit: microseconds
##    expr     min       lq      mean   median      uq      max neval cld
##  scalar 34428.2 36966.35 44244.423 38257.80 40968.6 407679.6   100   b
##  matrix    80.8    82.45   148.764   156.45   177.8    721.7   100  a
Метки: