Как векторизовать этот цикл? Умножьте две матрицы, сохраните информацию, сделайте это много раз без зацикливания

Предположим (маленькие числа в этом примере) у меня есть массив,

3 x 14 x 5

назови это

set.seed(1)
dfarray=array(rnorm(5*3*14,0,1),dim=c(3,14,5))

У меня есть матрица, которая соответствует этому и является

39 (which is 13*3) x 14

Назовите эту матрицу:

dfmat = matrix(rnorm(13*3*14,0,1),39,14)
dfmat = cbind(dfmat,rep(1:3,13))
dfmat = dfmat[order(dfmat [,15]),]
colnames(dfmat)[15]='unit'

Я хочу запустить этот цикл:

 costs = c(0.45, 2.11, 1.05, 1.44, 0.88, 2.30, 1.96, 1.76, 2.06, 1.54, 1.69,1.75,0)
    p = c(1,2,3,1,4,3,2,1,4,1,3,4,0)
    profit=numeric(0)
    for(i in 1:3){
            j=13
            beta = dfarray[i,,]
            Xt = dfmat [which(dfmat [,'unit']==i),1:14]    #this takes a set of 13, Xt is 13x14

            Xbeta = exp( Xt %*% beta )
            iota = c(rep(1, j))
            denom = iota%*%Xbeta
            Prob =  (Xbeta/ (iota%*%denom))
            Eprob = rowSums(Prob)/5  #the 5 coming from the last dim of array
            profit = c(profit,sum((p-costs)*Eprob))

        }


     sum(profit)

Я не могу придумать способ векторизовать ту часть, которую обходит цикл, вызывая

beta = dfarray[i,,]
Xt = dfmat [which(dfmat [,'unit']==i),]   #this takes a set of 13, Xt is 13x14

r loops apply matrix vectorization
person wolfsatthedoor    schedule 09.03.2015    source источник
comment
Можно ли было бы хранить dfmat также в трехмерном массиве?   -  person cryo111    schedule 09.03.2015
comment
Да, было бы, есть ли способ сделать умножение трехмерных массивов, которое будет работать? Я редактирую код, чтобы его можно было воспроизвести прямо сейчас, так что терпите меня!   -  person wolfsatthedoor    schedule 09.03.2015
comment
Взгляните на пакет tensorA. Он предоставляет функцию mul.tensor. Я предполагаю, что вам нужно использовать аргумент by в этой функции, чтобы перебрать ваш индекс i.   -  person cryo111    schedule 09.03.2015
comment
Я разберусь с этим прямо сейчас, теперь код можно воспроизводить, кстати!   -  person wolfsatthedoor    schedule 09.03.2015
comment
Функция выглядит как mul.tensor(X,i=c(),Y,j=i,by=NULL). X и Y будут вашими трехмерными массивами (в этом случае тензорными объектами). С помощью i и j вы можете определить индексы, которые используются для умножения матриц, а аргумент by - это ваша переменная цикла (i в вашем случае).   -  person cryo111    schedule 09.03.2015
comment
Векторизация для ускорения или удобочитаемости / обслуживания?   -  person Khashaa    schedule 09.03.2015
comment
Для ускорения. Цикл слишком медленный (у меня что-то вроде 20000 вместо 100)   -  person wolfsatthedoor    schedule 09.03.2015
comment
Тогда вы также можете использовать пакеты Rcpp и _2 _...   -  person cryo111    schedule 09.03.2015
comment
То есть переписать на C ++ правда?   -  person wolfsatthedoor    schedule 09.03.2015
comment
Mul.tensor не быстрее цикла?   -  person wolfsatthedoor    schedule 10.03.2015
comment
Написание кода RcppArmadillo должно быть довольно простым. Он имеет синтаксис, похожий на Matlab.   -  person Khashaa    schedule 10.03.2015
comment
Mul.tensor должно быть быстрее. Да, Rcpp использует C ++. Это займет некоторое время, но вместе с RcppArmadillo это действительно полезно, если вы делаете много вещей, которые нелегко векторизовать.   -  person cryo111    schedule 10.03.2015
comment
Собственно, я только что проверил mul.tensor - вроде не быстрее цикла. Не ожидал этого.   -  person cryo111    schedule 10.03.2015
comment
Казалось странным, что ваш воспроизводимый пример возвращает 104. Оказалось, что независимо от генерации данных sum(Eprob) равно 1, поэтому я немного озадачен тем, чего вы пытаетесь достичь.   -  person Khashaa    schedule 10.03.2015
comment
Я избавился от некоторых манипуляций с функциональными формами, но их легко добавить обратно. Я вынул для упрощения. Как вы думаете, мне стоит добавить еще раз?   -  person wolfsatthedoor    schedule 10.03.2015
comment
Есть ли какая-то особая причина хранить dfmat как именованную матрицу, а не список матриц?   -  person Khashaa    schedule 10.03.2015
comment
Основным узким местом является поиск имени Xt = dfmat [which(dfmat [,'unit']==i),] . Если вместо этого вы сохранитеdfmat как список матриц, а затем получите доступ к его элементам как Xt = dfmat [[i]], это будет в несколько раз быстрее.   -  person Khashaa    schedule 10.03.2015

Ответы (1)


arrow_upward
2
arrow_downward

Чтобы пояснить мои замечания в колонке комментариев, предположим, что у нас есть dfmat как список матриц. Почти всегда легче работать со списком матриц, чем с одной большой именованной матрицей. Также, если вы хотите полностью векторизовать приведенное здесь решение, вы можете получить блочную диагональную матрицу, используя bdiag из пакета Matrix, который действует со списками.

set.seed(1)
dfarray=array(rnorm(5*3*14,0,1),dim=c(3,14,5))
# dfmats as a list of matrices
dfmats <- lapply(1:3, function(i)matrix(rnorm(13*14), nrow=13))

Умножение на iota равно colSums или rowSums, поэтому мы можем упростить операцию, как в f.

f <- function(Xbeta) rowSums(Xbeta / matrix(colSums(Xbeta), nrow=nrow(Xbeta), ncol=ncol(Xbeta), byrow=T)) / ncol(Xbeta)
#profits is written as a function for benchmarking
#cost and p are ignored as they can be easily added back in.
profits <- function(){ 
    Xbetas <- lapply(seq_len(dim(dfarray)[1]), function(i) exp(dfmats[[i]] %*% dfarray[i,,]))
    Eprobs <- lapply(Xbetas, f)
    unlist(Eprobs)
}

И твой подход

profits1 <- function(){
    profit=numeric(0)
    for(i in 1:dim(dfarray)[1]){
        j=13
        beta = dfarray[i,,]
        Xt = dfmat [which(dfmat [,'unit']==i),1:14]    #this takes a set of 13, Xt is 13x14
        
        Xbeta = exp( Xt %*% beta )
        iota = c(rep(1, j))
        denom = iota%*%Xbeta
        deno <- colSums(Xbeta)
        s <- iota%*%denom
        Prob =  (Xbeta/ s)
        Eprob = rowSums(Prob)/dim(dfarray)[3]  #the 100 coming from the last dim of array
        profit = c(profit,Eprob)
        
    }
    return(profit)
}
dfmat <- do.call(rbind, dfmats)
dfmat <- cbind(dfmat,rep(1:3, each=13))
colnames(dfmat)[15]='unit'

Убедитесь, что они дают одинаковые результаты

all.equal(profits(), profits1())
[1] TRUE

Контрольный показатель

Я запустил это на экземпляре AWS EC2 free-tie, доступ к которому можно получить через http://www.louisaslett.com/RStudio_AMI/ < / а>.

dfarray=array(rnorm(100*10000*14,0,1),dim=c(10000,14,100))
dfmats <- lapply(1:10000, function(i)matrix(rnorm(13*14), nrow=13))

Из исходной конструкции вы можете преобразовать dfmat в список dfmats как dfmats <- lapply(1:3, function(i)dfmat[which(dfmat [,'unit']==i),1:14]), но это очень затратное преобразование. Создание dfmat из dfmats обходится дешевле.

dfmat <- do.call(rbind, dfmats)
dfmat <- cbind(dfmat,rep(1:10000, each=13))
colnames(dfmat)[15]='unit'

Обратите внимание на исключительное ускорение использования list и опасность ужасной стоимости поиска имени.

system.time(a1 <- profits1())
#   user  system elapsed 
#250.885   4.442 255.394 
system.time(a <- profits())
#   user  system elapsed 
#  2.717   0.429   3.167 
all.equal(a, a1)
#[1] TRUE

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

person Khashaa    schedule 15.03.2015