Я хотел бы ускорить выполнение функции R. Для каждого столбца в матрице 'A' найдите индекс (не сам), для которого его произведение с другим элементом вектора и соответствующим элементом в симметричной корреляционной матрице R является максимальным.
В настоящее время существует некоторая избыточность в вычислении внешнего произведения, поскольку оно без необходимости генерирует полную матрицу. Кроме того, в идеале цикл (то есть «применить») должен быть векторизован.
Пример данных ниже.
set.seed(123)
A <- matrix(rexp(30000, rate=.1), nrow=3000, ncol=2000)/100
R <- matrix( runif (10000), 3000 , 3000 )
diag(R) <- 1
R[upper.tri(R)] <- R[lower.tri(R)]
function_which_is_too_slow <- function(index){
aar <- outer(A[,index], A[,index]) * R
diag(aar) <- 0
return(max.col(aar, 'first'))
}
out <- apply(matrix(1:dim(A)[2]), 1, function_which_is_too_slow)
Спасибо, теперь воспроизводимо.
apply
вряд ли будет главной проблемой, если у вас всего 2000 колонок. Я не уверен, что вы могли бы добиться большего с базой R. Описание вашей проблемы, кажется, предполагает, что вам нужно будет посмотреть как минимум 2000 * 3000 * 2999/2 = 9 миллиардов продуктов. В настоящее время вы смотрите на продукты 2000 * 3000 ^ 2, что в 2 раза больше, но вы делаете это с помощью оптимизированных функций, таких как outer
и max.col
. Если время - проблема, это кажется естественным вариантом использования для rcpp
.
Вот ваш код как базовая линия с меньшим размером проблемы:
set.seed(123)
A <- matrix(rexp(30000, rate=.1), nrow=3000, ncol=40)/100
R <- matrix( runif (10000), 3000 , 3000 )
diag(R) <- 1
R[upper.tri(R)] <- R[lower.tri(R)]
function_which_is_too_slow <- function(index){
aar <- outer(A[,index], A[,index]) * R
diag(aar) <- 0
return(max.col(aar, 'first'))
}
system.time(apply(matrix(1:dim(A)[2]), 1, function_which_is_too_slow))
#> User System verstrichen
#> 12.001 11.388 10.348
Обнулять диагональ каждый раз не нужно, если мы используем копию корреляционной матрицы с нулевой диагональю. Использование lapply
вместо apply
выглядит лучше:
Rp <- R
diag(Rp) <- 0
faster_function <- function(index){
aar <- outer(A[,index], A[,index]) * Rp
return(max.col(aar, 'first'))
}
system.time(lapply(1:ncol(A), faster_function))
#> User System verstrichen
#> 11.156 10.306 8.334
Мы также можем использовать RcppArmadillo для выполнения тех же вычислений на C++.
Rcpp::cppFunction(code = "
arma::uvec arma_function(const arma::mat &A, const arma::mat &Rp, int index) {
arma::mat aar = A.col(index-1) * A.col(index-1).t() % Rp;
return index_max(aar, 1) + 1;
}
", depends = "RcppArmadillo")
system.time(lapply(1:ncol(A), arma_function, A = A, Rp = Rp))
#> User System verstrichen
#> 7.382 10.578 4.874
И мы можем распараллелить вычисления, хотя RcppArmadillo уже использует OpenMP, если он доступен:
system.time(parallel::mclapply(1:ncol(A), arma_function, A = A, Rp = Rp))
#> User System verstrichen
#> 0.008 0.010 3.903
В целом, примерно в 3 раза быстрее, что немного.
Я знаю, что это не играет роли, но для хорошей практики используйте
set.seed
при использовании таких функций, какrunif
и т. д.