Быстро умножить все комбинации матриц в списке - PullRequest
Новая
       6

Быстро умножить все комбинации матриц в списке

0 голосов
/ 19 февраля 2019

У меня есть симметричная матрица Siginv и список Z, содержащий N матрицы размера TxK.

У меня есть три подхода для вычисления того, что мне нужно (см. Ниже).Первый буквально принимает формулу и работает медленно, как и ожидалось.Есть ли более быстрые способы делать то, что я делаю? 

library(microbenchmark)

K <- 2
N <- 50
Siginv <- matrix(rnorm(N^2), ncol=N)
Siginv[lower.tri(Siginv)] <- t(Siginv)[lower.tri(Siginv)] # just some symmetric matrix

Tdim <- 400
Z <- replicate(N, matrix(rnorm(Tdim*K), ncol=K), simplify = F)

microbenchmark({
  I <- diag(Tdim)
  Z.m <- do.call(rbind, Z)
  meat.mat.GLS.kp <- t(Z.m)%*%(Siginv%x%I)%*%Z.m
}, {
  combs <- expand.grid(1:N, 1:N)
  cprods.GLS <- mapply(function(i,j) Siginv[j,i]*t(Z[[i]])%*%Z[[j]], combs[,1], combs[,2], SIMPLIFY = F)
  meat.mat.GLS <- Reduce("+", cprods.GLS)
}, {
  combs <- expand.grid(1:N, 1:N)
  cprods.GLS2 <- mapply(function(i,j) Siginv[j,i]*crossprod(Z[[i]],Z[[j]]), combs[,1], combs[,2], SIMPLIFY = F)
  meat.mat.GLS2 <- Reduce("+", cprods.GLS2)
}, times=5)

all.equal(meat.mat.GLS.kp, meat.mat.GLS, meat.mat.GLS2) # TRUE

В настоящее время подходы сравниваются следующим образом: 

        min         lq       mean     median         uq        max neval cld
 4499.66564 4911.42674 4874.35170 4958.81553 4977.55873 5024.29187     5   b
   23.03861   23.09293   23.82407   23.29574   24.04696   25.64611     5  a 
   12.92261   13.08275   13.54088   13.15898   13.80212   14.73794     5  a
...