Ускорить умножение векторов

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

Я знаю, что это не играет роли, но для хорошей практики используйте set.seed при использовании таких функций, как runif и т. д.

Sotos 10.08.2018 11:11

Спасибо, теперь воспроизводимо.

user7105520 10.08.2018 11:16

apply вряд ли будет главной проблемой, если у вас всего 2000 колонок. Я не уверен, что вы могли бы добиться большего с базой R. Описание вашей проблемы, кажется, предполагает, что вам нужно будет посмотреть как минимум 2000 * 3000 * 2999/2 = 9 миллиардов продуктов. В настоящее время вы смотрите на продукты 2000 * 3000 ^ 2, что в 2 раза больше, но вы делаете это с помощью оптимизированных функций, таких как outer и max.col. Если время - проблема, это кажется естественным вариантом использования для rcpp.

John Coleman 10.08.2018 15:09
Стоит ли изучать PHP в 2023-2024 годах?
Стоит ли изучать PHP в 2023-2024 годах?
Привет всем, сегодня я хочу высказать свои соображения по поводу вопроса, который я уже много раз получал в своем сообществе: "Стоит ли изучать PHP в...
Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией
Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией
В JavaScript одним из самых запутанных понятий является поведение ключевого слова "this" в стрелочной и обычной функциях.
Приемы CSS-макетирования - floats и Flexbox
Приемы CSS-макетирования - floats и Flexbox
Здравствуйте, друзья-студенты! Готовы совершенствовать свои навыки веб-дизайна? Сегодня в нашем путешествии мы рассмотрим приемы CSS-верстки - в...
Тестирование функциональных ngrx-эффектов в Angular 16 с помощью Jest
В системе управления состояниями ngrx, совместимой с Angular 16, появились функциональные эффекты. Это здорово и делает код определенно легче для...
Концепция локализации и ее применение в приложениях React ⚡️
Концепция локализации и ее применение в приложениях React ⚡️
Локализация - это процесс адаптации приложения к различным языкам и культурным требованиям. Это позволяет пользователям получить опыт, соответствующий...
Пользовательский скаляр GraphQL
Пользовательский скаляр GraphQL
Листовые узлы системы типов GraphQL называются скалярами. Достигнув скалярного типа, невозможно спуститься дальше по иерархии типов. Скалярный тип...
1
3
154
1

Ответы 1

Вот ваш код как базовая линия с меньшим размером проблемы:

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 раза быстрее, что немного.

Другие вопросы по теме