Умножение многих матриц в R - PullRequest
Новая
       53

Умножение многих матриц в R

0 голосов
/ 13 ноября 2018

Я хочу умножить несколько матриц одинакового размера на начальный вектор.В приведенном ниже примере p.state является вектором m элементов, а tran.mat является списком, где каждый элемент является m x m матрицей.

for (i in 1:length(tran.mat)){
  p.state <- p.state %*% tran.mat[[i]]
}

Приведенный выше код дает правильный ответ, но может быть медленным, если length(tran.mat) велико.Мне было интересно, есть ли более эффективный способ сделать это?

Ниже приведен пример с m=3 и length(mat)=10, который может генерировать это: 

p.state <- c(1,0,0)
tran.mat<-lapply(1:10,function(y){apply(matrix(runif(9),3,3),1,function(x){x/sum(x)})})

for (i in 1:length(tran.mat)){
  p.state <- p.state %*% tran.mat[[i]]
}

print(p.state)

NB: tran.mat не обязательно должен быть списком, он просто записан как единое целое.

Редактировать после нескольких комментариев:

Reduce полезно, когда m мало.Однако при m=6 выход из цикла выполнил оба вышеупомянутых решения.библиотека (rbenchmark) 

p.state1 <- p.state <- c(1,0,0,0,0,0)
tran.mat<-lapply(1:10000,function(y){t(apply(matrix(runif(36),6,6),1,function(x){x/sum(x)}))})

tst<-do.call(c, list(list(p.state), tran.mat))

benchmark(
  'loop' = {
    for (i in 1:length(tran.mat)){
      p.state <- p.state %*% tran.mat[[i]]
    }
  },
  'reduce' = {
    p.state1 %*% Reduce('%*%', tran.mat)
   },
  'reorder' = {
    Reduce(`%*%`,tran.mat,p.state1)
  }

)

В результате 

        test replications elapsed relative user.self sys.self user.child     sys.child
   1    loop          100    0.87    1.000      0.87        0         NA        NA
   2  reduce          100    1.41    1.621      1.39        0         NA        NA
   3 reorder          100    1.00    1.149      1.00        0         NA        NA

Ответы [ 3 ]

0 голосов
/ 13 ноября 2018

Я не уверен, что это будет быстрее, но короче: 

prod <- Reduce("%*%", L)

all.equal(prod, L[[1]] %*% L[[2]] %*% L[[3]] %*% L[[4]])
## [1] TRUE

Примечание

Мы использовали этот тестовый вход: 

m <- matrix(1:9, 3)
L <- list(m^0, m, m^2, m^3)
0 голосов
/ 17 ноября 2018

Я собираюсь использовать функцию из пакета Rfast , чтобы сократить время выполнения умножения. К сожалению, время цикла не может быть сокращено.

Функция с именем Rfast::eachcol.apply - отличное решение для ваших целей. Ваше умножение также является функцией crossprod, но для нашей цели оно медленное.

Вот некоторые вспомогательные функции: 

mult.list<-function(x,y){
    for (xm in x){
        y <- y %*% xm
    }
    y
}

mult.list2<-function(x,y){
    for (xm in x){
        y <- Rfast::eachcol.apply(xm,y,oper="*",apply="sum")
    }
    y
}

Вот пример: 

x<-list()
y<-rnomr(1000)
for(i in 1:100){
    x[[i]]<-Rfast::matrnorm(1000,1000)
}


microbenchmark::microbenchmark(R=a<-mult.list(x,y),Rfast=b<-mult.list2(x,y),times = 10)
 Unit: milliseconds
     expr        min         lq        mean     median         uq        max neval
        R 410.067525 532.176979 633.3700627 649.155826 699.721086 916.542414    10
    Rfast 239.987159 251.266488 352.1951486 276.382339 458.089342 741.340268    10

all.equal(as.numeric(a),as.numeric(b))
[1] TRUE

Аргумент oper предназначен для операции над каждым элементом, а применяет для операции над каждым столбцом. В больших матрицах должно быть быстро. Я не мог проверить это на своем ноутбуке для больших матриц.

0 голосов
/ 13 ноября 2018

Более быстрый способ - использовать Reduce() для последовательного умножения матриц в списке матриц.

Таким образом вы можете получить примерно 4-кратное ускорение. Ниже приведен пример протестированного кода с 1000 элементами в списке вместо 10, чтобы упростить процесс улучшения производительности.

код 

library(rbenchmark)

p.state <- c(1,0,0)
tran.mat<-lapply(1:1000,function(y){apply(matrix(runif(9),3,3),1,function(x){x/sum(x)})})

benchmark(
  'loop' = {
    for (i in 1:length(tran.mat)){
      p.state <- p.state %*% tran.mat[[i]]
    }
  },
  'reduce' = {
    p.state %*% Reduce('%*%', tran.mat)
  }
)

выход 

    test replications elapsed relative user.self sys.self user.child sys.child
1   loop          100    0.23    3.833      0.23        0         NA        NA
2 reduce          100    0.06    1.000      0.07        0         NA        NA

Вы можете видеть, что метод сокращения примерно в 3,8 раза быстрее.

...