Можно ли векторизовать расчет, который зависит от предыдущих элементов?

Ответ 1

ifelse векторизован и существует немного штрафа, если вы используете его по одному элементу за раз в цикле for. В вашем примере вы можете получить довольно хорошее ускорение, используя if вместо ifelse.

fun1 <- function(z) {
  for(i in 2:NROW(z)) {
    z[i] <- ifelse(z[i-1]==1, 1, 0)
  }
  z
}

fun2 <- function(z) {
  for(i in 2:NROW(z)) {
    z[i] <- if(z[i-1]==1) 1 else 0
  }
  z
}

z <- c(1,1,0,0,0,0)
identical(fun1(z),fun2(z))
# [1] TRUE
system.time(replicate(10000, fun1(z)))
#   user  system elapsed 
#   1.13    0.00    1.32
system.time(replicate(10000, fun2(z)))
#   user  system elapsed 
#   0.27    0.00    0.26 

Вы можете получить дополнительную прибавку скорости из fun2, скомпилировав ее.

library(compiler)
cfun2 <- cmpfun(fun2)
system.time(replicate(10000, cfun2(z)))
#   user  system elapsed 
#   0.11    0.00    0.11

Итак, 10-кратное ускорение без векторизации. Как говорили другие (и некоторые из них проиллюстрированы) есть способы для векторизации вашего примера, но это может не привести к вашей реальной проблеме. Надеюсь, это достаточно общее, чтобы быть применимым.

Функция filter может быть полезна и вам, если вы можете понять, как выразить свою проблему с точки зрения процесса авторегрессии или скользящего среднего.

Ответ 2

Это хороший и простой пример, где Rcpp может сиять.

Итак, сначала переработаем функции 1 и 2 и их скомпилированные копии:

library(inline)
library(rbenchmark)
library(compiler)

fun1 <- function(z) {
  for(i in 2:NROW(z)) {
    z[i] <- ifelse(z[i-1]==1, 1, 0)
  }
  z
}
fun1c <- cmpfun(fun1)


fun2 <- function(z) {
  for(i in 2:NROW(z)) {
    z[i] <- if(z[i-1]==1) 1 else 0
  }
  z
}
fun2c <- cmpfun(fun2)

Мы очень легко пишем вариант Rcpp:

funRcpp <- cxxfunction(signature(zs="numeric"), plugin="Rcpp", body="
  Rcpp::NumericVector z = Rcpp::NumericVector(zs);
  int n = z.size();
  for (int i=1; i<n; i++) {
    z[i] = (z[i-1]==1.0 ? 1.0 : 0.0);
  }
  return(z);
")

Здесь используется пакет inline для компиляции, загрузки и привязки пятистрочного интерфейса на лету.

Теперь мы можем определить нашу дату теста, которую мы делаем немного дольше, чем оригинал (так как только запуск оригинала слишком мало раз приводит к неизмеримым временам):

R> z <- rep(c(1,1,0,0,0,0), 100)
R> identical(fun1(z),fun2(z),fun1c(z),fun2c(z),funRcpp(z))
[1] TRUE
R> 

Все ответы рассматриваются как идентичные.

Наконец, мы можем сравнить:

R> res <- benchmark(fun1(z), fun2(z),
+                  fun1c(z), fun2c(z),
+                  funRcpp(z),
+                  columns=c("test", "replications", "elapsed", 
+                            "relative", "user.self", "sys.self"),
+                  order="relative",
+                  replications=1000)
R> print(res)
        test replications elapsed relative user.self sys.self
5 funRcpp(z)         1000   0.005      1.0      0.01        0
4   fun2c(z)         1000   0.466     93.2      0.46        0
2    fun2(z)         1000   1.918    383.6      1.92        0
3   fun1c(z)         1000  10.865   2173.0     10.86        0
1    fun1(z)         1000  12.480   2496.0     12.47        0

Скомпилированная версия выигрывает почти в 400 раз против лучшей версии R и почти 100 против ее байтового скомпилированного варианта. Для функции 1 байт-компиляция имеет гораздо меньшее значение, и оба варианта отслеживают С++ в два раза более чем на две тысячи.

Потребовалось около одной минуты, чтобы написать версию на С++. Ускорение скорости показывает, что это была потраченная минуту.

Для сравнения, вот результат для исходного короткого вектора, который называется чаще:

R> z <- c(1,1,0,0,0,0)
R> res2 <- benchmark(fun1(z), fun2(z),
+                  fun1c(z), fun2c(z),
+                  funRcpp(z),
+                  columns=c("test", "replications", 
+                            "elapsed", "relative", "user.self", "sys.self"),
+                  order="relative",
+                  replications=10000)
R> print(res2)
        test replications elapsed  relative user.self sys.self
5 funRcpp(z)        10000   0.046  1.000000      0.04        0
4   fun2c(z)        10000   0.132  2.869565      0.13        0
2    fun2(z)        10000   0.271  5.891304      0.27        0
3   fun1c(z)        10000   1.045 22.717391      1.05        0
1    fun1(z)        10000   1.202 26.130435      1.20        0

Качественное ранжирование не изменилось: доминирует версия Rcpp, функция 2 - вторая. с байтовой скомпилированной версией примерно в два раза быстрее, чем простой вариант R, но все же почти в три раза медленнее, чем версия С++. И относительная разница ниже: относительно говоря, служебная информация вызова функции меньше, а фактическая петля имеет значение больше: С++ получает большее преимущество от фактических операций цикла в более длинных векторах. То, что это важный результат, поскольку он предполагает, что больше данных реального размера, скомпилированная версия может получить большую выгоду.

Отредактировано для исправления двух небольших недочетов в примерах кода. И снова отредактирован с благодарностью Джошу, чтобы поймать ошибку установки относительно fun2c.

Ответ 3

Я думаю, что это обман, а не обобщаемый, но: согласно правилам, которые вы указали выше, любое появление 1 в векторе сделает все последующие элементы 1 (по рекурсии: z[i] равно 1, если z[i-1] равно 1, поэтому z[i] будет установлено в 1, если z[i-2] равно 1 и т.д.). В зависимости от того, что вы действительно хотите сделать, может быть такое рекурсивное решение, если вы тщательно об этом подумаете...

z <- c(1,1,0,0,0,0)
first1 <- min(which(z==1))
z[seq_along(z)>first1] <- 1

edit: это неправильно, но я оставляю его, чтобы признать свои ошибки. Основываясь на немного игры (и меньше думая), я думаю, что фактическое решение этой рекурсии более симметрично и даже проще:

rep(z[1],length(z))

Тестовые случаи:

z <- c(1,1,0,0,0,0)
z <- c(0,1,1,0,0,0)
z <- c(0,0,1,0,0,0)

Ответ 4

Проверьте rollapply функцию zoo.

Я не очень хорошо знаком с этим, но я думаю делает то, что вам нужно:

> c( 1, rollapply(z,2,function(x) x[1]) )
[1] 1 1 1 1 1 1

Я как бы блокирую его, используя окно из 2, а затем только используя первый элемент этого окна.

Для более сложных примеров вы можете выполнить некоторые вычисления на x [1] и вместо этого вернуть это.

Ответ 5

Иногда вам просто нужно думать об этом совершенно по-другому. То, что вы делаете, создает вектор, в котором каждый элемент будет таким же, как первый, если он равен 1 или 0. В противном случае.

z <- c(1,1,0,0,0,0)
if (z[1] != 1) z[1] <- 0
z[2:length(z)] <- z[1]

Ответ 6

Это также можно сделать с помощью "apply" с использованием исходного вектора и запаздывающей версии вектора в качестве составных столбцов кадра данных.