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

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

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 ответа


Более быстрый способ - использовать 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 раза быстрее.


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

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)

Я собираюсь использовать функцию из пакета 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

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


Есть идеи?

10000