Векторизация в R сложного алгоритма

#r #performance #vectorization

#r #Производительность #векторизация

Вопрос:

Мой код вызывает приведенную ниже функцию много раз. Я использовал Rprof, чтобы выяснить, что он занимает 1/3 времени выполнения.

Я слышал, что R можно ускорить с помощью векторизации. Однако алгоритм, который я использую, обращается к двум индексам в двух разных списках, поэтому я не вижу, как будет работать какая-либо из функций *ply .

Есть ли что-нибудь еще, что можно было бы сделать для оптимизации этого?

n — целое число, u и v — списки комплексных чисел

 psi <- function(n, u, v)
{
    psi = complex(real = 0, imaginary = 0)

    for (i in 1 : (n - 1))
    {
        for (j in (i   1) : n)
        {
            psi = psi   log(u[i] * v[j] - u[j] * v[i])
        }
    }

    return (psi * 3)
}
  

Ответ №1:

Я считаю, что следующее может решить вашу проблему. Я выбрал этот метод выражения, потому что считаю его очень прозрачным, учитывая то, что вы уже сделали. По сути, я только что заранее сгенерировал вектор i и j и использовал их так же, как и вы, но одновременно как векторы (скопировал и вставил ваше уравнение).

 psi <- function(n,u,v) {
    j <- combn(n, 2) #I'll overwrite j later to conserve memory
    i <- j[1,]
    j <- j[2,]
    sum(log(u[i] * v[j] - u[j] * v[i])) * 3
}
  

Это будет намного быстрее, чем ваши циклы, если n не очень большое число, потому что оно должно генерировать и поддерживать в памяти весь набор комбинаций, который вы генерируете с помощью своих циклов одновременно. Но в тот момент, когда это приведет к сбою, вам в любом случае придется долго вычислять.

Чтобы показать эквивалентность, рассмотрим результаты сравнения ваших методов генерации ваших индексов и настоящего.

 n <- 5
combn(n, 2)

for (i in 1 : (n - 1))
{
    for (j in (i   1) : n)
    {
        print(c(i, j))
    }
}
  

Ответ №2:

Если я правильно прочитал ваш код, для случая n=4 , когда ваше суммирование сводится к

 log(u[1]*v[2] - u[2]*v[1])  
log(u[1]*v[3] - u[3]*v[1])  
log(u[1]*v[4] - u[4]*v[1])  

log(u[2]*v[3] - u[3]*v[2])  
log(u[2]*v[4] - u[4]*v[2])  

log(u[3]*v[4] - u[4]*v[3])
  

на самом деле это просто все возможные комбинации из n, выберите 2. Вы можете создать такой список индексов с combn() помощью .

 psi <- function(n,u,v) {
    mx <- function(p) {a<-p[1]; b<-p[2]; log(u[a]*v[b]-u[b]*v[a])}
    sum(combn(n,2, FUN=mx)) * 3
}
  

Здесь внутренняя mx функция выполняет вычисления для каждой пары индексов. Кажется, это заметно быстрее. Я думаю, вы хотите убедиться, что они получают одинаковые результаты. Мне было трудно генерировать тестовые входные данные, которые не вызывали log гнева.