Действительно ли семейство "* apply" не векторизовано?

Ответ 1

Прежде всего, в вашем примере вы делаете тесты на "data.frame", что не справедливо для colMeans, apply и "[.data.frame", поскольку у них есть служебные данные:

system.time(as.matrix(m))  #called by `colMeans` and `apply`
#   user  system elapsed 
#   1.03    0.00    1.05
system.time(for(i in 1:ncol(m)) m[, i])  #in the `for` loop
#   user  system elapsed 
#  12.93    0.01   13.07

На матрице изображение немного отличается:

mm = as.matrix(m)
system.time(colMeans(mm))
#   user  system elapsed 
#   0.01    0.00    0.01 
system.time(apply(mm, 2, mean))
#   user  system elapsed 
#   1.48    0.03    1.53 
system.time(for(i in 1:ncol(mm)) mean(mm[, i]))
#   user  system elapsed 
#   1.22    0.00    1.21

Заменяя основную часть вопроса, основное различие между lapply/mapply/etc и прямыми R-петлями - это то, где выполняется цикл. Как отмечает Роланд, как циклы C, так и R должны оценивать функцию R на каждой итерации, которая является наиболее дорогостоящей. Очень быстрые функции C - это те, которые делают все на C, поэтому, я думаю, это должно быть то, о чем "векторизован"?

Пример, где мы находим среднее значение в каждом из элементов "списка":

( EDIT 11 мая '16. Я считаю, что пример с поиском "среднего" не является хорошей настройкой различий между оценкой функции R итеративно и скомпилированным кодом, (1) из-за особенность R-алгоритма на "числовых" по сравнению с простыми sum(x) / length(x) и (2) должна иметь смысл тестировать "список" с помощью length(x) >> lengths(x). Таким образом, "средний" пример переносится в конец и заменить другим.)

В качестве простого примера мы могли бы рассмотреть нахождение противоположности каждого элемента length == 1 в "списке":

В файле tmp.c:

#include <R.h>
#define USE_RINTERNALS 
#include <Rinternals.h>
#include <Rdefines.h>

/* call a C function inside another */
double oppC(double x) { return(ISNAN(x) ? NA_REAL : -x); }
SEXP sapply_oppC(SEXP x)
{
    SEXP ans = PROTECT(allocVector(REALSXP, LENGTH(x)));
    for(int i = 0; i < LENGTH(x); i++) 
        REAL(ans)[i] = oppC(REAL(VECTOR_ELT(x, i))[0]);

    UNPROTECT(1);
    return(ans);
}

/* call an R function inside a C function;
 * will be used with 'f' as a closure and as a builtin */    
SEXP sapply_oppR(SEXP x, SEXP f)
{
    SEXP call = PROTECT(allocVector(LANGSXP, 2));
    SETCAR(call, install(CHAR(STRING_ELT(f, 0))));

    SEXP ans = PROTECT(allocVector(REALSXP, LENGTH(x)));     
    for(int i = 0; i < LENGTH(x); i++) { 
        SETCADR(call, VECTOR_ELT(x, i));
        REAL(ans)[i] = REAL(eval(call, R_GlobalEnv))[0];
    }

    UNPROTECT(2);
    return(ans);
}

И со стороны R:

system("R CMD SHLIB /home/~/tmp.c")
dyn.load("/home/~/tmp.so")

с данными:

set.seed(007)
myls = rep_len(as.list(c(NA, runif(3))), 1e7)

#a closure wrapper of `-`
oppR = function(x) -x

for_oppR = compiler::cmpfun(function(x, f)
{
    f = match.fun(f)  
    ans = numeric(length(x))
    for(i in seq_along(x)) ans[[i]] = f(x[[i]])
    return(ans)
})

Бенчмаркинг:

#call a C function iteratively
system.time({ sapplyC =  .Call("sapply_oppC", myls) }) 
#   user  system elapsed 
#  0.048   0.000   0.047 

#evaluate an R closure iteratively
system.time({ sapplyRC =  .Call("sapply_oppR", myls, "oppR") }) 
#   user  system elapsed 
#  3.348   0.000   3.358 

#evaluate an R builtin iteratively
system.time({ sapplyRCprim =  .Call("sapply_oppR", myls, "-") }) 
#   user  system elapsed 
#  0.652   0.000   0.653 

#loop with a R closure
system.time({ forR = for_oppR(myls, "oppR") })
#   user  system elapsed 
#  4.396   0.000   4.409 

#loop with an R builtin
system.time({ forRprim = for_oppR(myls, "-") })
#   user  system elapsed 
#  1.908   0.000   1.913 

#for reference and testing 
system.time({ sapplyR = unlist(lapply(myls, oppR)) })
#   user  system elapsed 
#  7.080   0.068   7.170 
system.time({ sapplyRprim = unlist(lapply(myls, `-`)) }) 
#   user  system elapsed 
#  3.524   0.064   3.598 

all.equal(sapplyR, sapplyRprim)
#[1] TRUE 
all.equal(sapplyR, sapplyC)
#[1] TRUE
all.equal(sapplyR, sapplyRC)
#[1] TRUE
all.equal(sapplyR, sapplyRCprim)
#[1] TRUE
all.equal(sapplyR, forR)
#[1] TRUE
all.equal(sapplyR, forRprim)
#[1] TRUE

(выполняется исходный пример среднего поиска):

#all computations in C
all_C = inline::cfunction(sig = c(R_ls = "list"), body = '
    SEXP tmp, ans;
    PROTECT(ans = allocVector(REALSXP, LENGTH(R_ls)));

    double *ptmp, *pans = REAL(ans);

    for(int i = 0; i < LENGTH(R_ls); i++) {
        pans[i] = 0.0;

        PROTECT(tmp = coerceVector(VECTOR_ELT(R_ls, i), REALSXP));
        ptmp = REAL(tmp);

        for(int j = 0; j < LENGTH(tmp); j++) pans[i] += ptmp[j];

        pans[i] /= LENGTH(tmp);

        UNPROTECT(1);
    }

    UNPROTECT(1);
    return(ans);
')

#a very simple `lapply(x, mean)`
C_and_R = inline::cfunction(sig = c(R_ls = "list"), body = '
    SEXP call, ans, ret;

    PROTECT(call = allocList(2));
    SET_TYPEOF(call, LANGSXP);
    SETCAR(call, install("mean"));

    PROTECT(ans = allocVector(VECSXP, LENGTH(R_ls)));
    PROTECT(ret = allocVector(REALSXP, LENGTH(ans)));

    for(int i = 0; i < LENGTH(R_ls); i++) {
        SETCADR(call, VECTOR_ELT(R_ls, i));
        SET_VECTOR_ELT(ans, i, eval(call, R_GlobalEnv));
    }

    double *pret = REAL(ret);
    for(int i = 0; i < LENGTH(ans); i++) pret[i] = REAL(VECTOR_ELT(ans, i))[0];

    UNPROTECT(3);
    return(ret);
')                    

R_lapply = function(x) unlist(lapply(x, mean))                       

R_loop = function(x) 
{
    ans = numeric(length(x))
    for(i in seq_along(x)) ans[i] = mean(x[[i]])
    return(ans)
} 

R_loopcmp = compiler::cmpfun(R_loop)


set.seed(007); myls = replicate(1e4, runif(1e3), simplify = FALSE)
all.equal(all_C(myls), C_and_R(myls))
#[1] TRUE
all.equal(all_C(myls), R_lapply(myls))
#[1] TRUE
all.equal(all_C(myls), R_loop(myls))
#[1] TRUE
all.equal(all_C(myls), R_loopcmp(myls))
#[1] TRUE

microbenchmark::microbenchmark(all_C(myls), 
                               C_and_R(myls), 
                               R_lapply(myls), 
                               R_loop(myls), 
                               R_loopcmp(myls), 
                               times = 15)
#Unit: milliseconds
#            expr       min        lq    median        uq      max neval
#     all_C(myls)  37.29183  38.19107  38.69359  39.58083  41.3861    15
#   C_and_R(myls) 117.21457 123.22044 124.58148 130.85513 169.6822    15
#  R_lapply(myls)  98.48009 103.80717 106.55519 109.54890 116.3150    15
#    R_loop(myls) 122.40367 130.85061 132.61378 138.53664 178.5128    15
# R_loopcmp(myls) 105.63228 111.38340 112.16781 115.68909 128.1976    15

Ответ 2

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

Цель векторной функции - исключить бухгалтерский учет, связанный с циклом for. Например, вместо:

means <- numeric(length(mtcars))
for (i in seq_along(mtcars)) {
  means[i] <- mean(mtcars[[i]])
}
sds <- numeric(length(mtcars))
for (i in seq_along(mtcars)) {
  sds[i] <- sd(mtcars[[i]])
}

Вы можете написать:

means <- vapply(mtcars, mean, numeric(1))
sds   <- vapply(mtcars, sd, numeric(1))

Это упрощает просмотр того, что же (входные данные) и что другое (функция, которую вы применяете).

Второе преимущество векторизации заключается в том, что for-loop часто записывается на C, а не в R. Это имеет значительные преимущества в производительности, но я не считаю это ключевым свойством векторизации. Векторизация в основном касается сохранения вашего мозга, а не сохранения работы компьютера.

Ответ 3

Я согласен с представлением Патрика Бернса о том, что это скорее loop hiding, а не кодирование векторов. Вот почему:

Рассмотрим этот фрагмент кода C:

  for (int я = 0; я < n; я ++) c [i] = a [i] + b [i]
Код>

Что мы хотели бы сделать, это совершенно ясно. Но , как выполнялась задача или как ее можно было выполнить, на самом деле не так. A for-loop по умолчанию является последовательной конструкцией. Он не сообщает, что и как можно сделать параллельно.

Наиболее очевидным способом является то, что код выполняется в последовательном порядке. Загрузите a [i] и b [i] для регистрации, добавления их, сохраните результат в c [i] и сделайте это для каждого i.

Однако современные процессоры имеют вектор или SIMD, который может работать с вектором данных во время той же инструкции при выполнении той же операции (например, добавление двух векторов, как показано выше). В зависимости от процессора/архитектуры может быть возможно добавить, например, четыре номера из a и b в соответствии с одной и той же инструкцией вместо одной за раз.

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

Обратите внимание, что это отличается от параллелизма кода, когда одновременно выполняются множественные вычисления.

Было бы здорово, если компилятор идентифицирует такие блоки кода, а автоматически их вектурирует, что является трудной задачей. Автоматическая векторизация кода является сложной темой исследований в области компьютерных наук. Но со временем компиляторы стали лучше. Вы можете проверить возможности автоматического векторизации GNU-gcc здесь. Аналогично для LLVM-clang здесь. И вы также можете найти некоторые тесты в последней ссылке по сравнению с gcc и ICC (компилятор Intel С++).

gcc (я нахожусь на v4.9)), например, не индексирует код автоматически при оптимизации уровня -O2. Поэтому, если бы мы выполнили код, показанный выше, он будет запускаться последовательно. Вот время для добавления двух целых векторов длиной 500 миллионов.

Нам нужно добавить флаг -ftree-vectorize или изменить оптимизацию до уровня -O3. (Обратите внимание, что -O3 выполняет другие дополнительные оптимизации). Флаг -fopt-info-vec полезен, поскольку он информирует, когда цикл был успешно проиндексирован).

  # компиляция с -O2, -ftree-vectorize и -fopt-info-vec
# test.c: 32: 5: примечание: цикл векторизован
# test.c: 32: 5: примечание: цикл с версией для векторизации из-за возможного сглаживания
# test.c: 32: 5: примечание: очистка петли для векторизации для улучшения выравнивания
Код>

Это говорит нам, что функция векторизована. Ниже приведены временные значения, сравнивающие как не-векторизованные, так и векторизованные версии на целых векторах длиной 500 миллионов:

  x = sample (100L, 500e6L, TRUE)
y = образец (100L, 500e6L, TRUE)
z = вектор ( "целое число", 500e6L) # вектор результата

# не-векторный, -O2
system.time(.Call( "Csum", x, y, z))
Прошла 1 пользовательская система
# 1,830 0,009 1,852

# vectorized с использованием флагов, показанных выше, на -O2
system.time(.Call( "Csum", x, y, z))
Прошла 1 пользовательская система
# 0,361 0,001 0,322

# оба результата проверяются на идентичность, возвращает TRUE
Код>

Эта часть может быть безопасно пропущена без потери непрерывности.

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

  #pragma omp simd
для (i = 0; я < n; я ++) c [i] = a [i] + b [i]
Код>

Таким образом, мы специально просим компилятор проинвестизировать его независимо от того, что. Нам нужно активировать расширения OpenMP, используя флаг времени компиляции -fopenmpкод>. Делая это:

  # timing с -O2 + OpenMP с simd
x = образец (100L, 500e6L, TRUE)
y = образец (100L, 500e6L, TRUE)
z = вектор ( "целое число", 500e6L) # вектор результата
system.time(.Call( "Cvecsum", x, y, z))
Прошла 1 пользовательская система
# 0,360 0,001 0,360
Код>

Это здорово! Это было протестировано с помощью gcc v6.2.0 и llvm clang v3.9.0 (оба установлены через homebrew, MacOS 10.12.3), оба из которых поддерживают OpenMP 4.0.

В этом смысле, хотя страница Википедии о программировании Array упоминает, что языки, которые работают на всей массивы обычно называют это векторизованными операциями, на самом деле это loop hiding IMO (если это фактически не векторизация).

В случае R код rowSums() или colSums() в C не использует идентификацию вектора IIUC; это всего лишь петля в C. То же самое касается lapply(). В случае apply(), он находится в R. Таким образом, все они скрываются в цикле.

Короче говоря, обертывание функции R:

просто написать for-loop в C!= векторный код вашего кода.
просто записывая for-loop в R!= векторный код вашего кода.

Intel Math Kernel Library (MKL), например, реализует векторизованные формы функций.

НТН

Ссылки:

• Обсуждение Джеймсом Рейндерсом, Intel (этот ответ в основном является попыткой обобщить это отличный разговор)

Ответ 4

Поэтому, чтобы суммировать большие ответы/комментарии до некоторого общего ответа и предоставить некоторый фон: R имеет 4 типа петель (из не-векторизованного в векторный порядок)

1. R for цикла, который неоднократно вызывает R-функции в каждой итерации (не-векторизован)
2. C, который неоднократно вызывает R-функции в каждой итерации (без векторизации)
3. C, который вызывает функцию R только один раз (несколько векторизованный)
4. Простой C-цикл, который вообще не вызывает никакой функции R и использует его собственные скомпилированные функции (в векторном виде)

Таким образом, семейство *apply является вторым типом. Кроме apply который больше первого типа

Вы можете понять это из комментария в своем исходном коде

/*.Internal(lapply (X, FUN)) */

/* Это особый. Внутри, поэтому есть неоцененные аргументы. это
вызывается из закрывающей оболочки, поэтому X и FUN являются обещаниями. FUN должно быть не оценено для использования, например, в bquote. */

Это означает, что lapply код C принимает неоцененную функцию из R, а затем оценивает ее внутри самого кода C. Это в основном разница между lapply .Internal вызов

.Internal(lapply(X, FUN))

У которого есть аргумент FUN который содержит функцию R

И colMeans .Internal вызов, у которого нет аргумента FUN

.Internal(colMeans(Re(x), n, prod(dn), na.rm))

colMeans, в отличие от lapply точно знает, какую функцию он должен использовать, поэтому он вычисляет среднее внутри внутри кода C.

Вы можете четко видеть процесс оценки функции R на каждой итерации в lapply C-кода

 for(R_xlen_t i = 0; i < n; i++) {
      if (realIndx) REAL(ind)[0] = (double)(i + 1);
      else INTEGER(ind)[0] = (int)(i + 1);
      tmp = eval(R_fcall, rho);   // <----------------------------- here it is
      if (MAYBE_REFERENCED(tmp)) tmp = lazy_duplicate(tmp);
      SET_VECTOR_ELT(ans, i, tmp);
   }

Чтобы подвести итог, lapply не векторизован, хотя он имеет два возможных преимущества перед lapply R for цикла

1. Доступ и назначение в цикле, по-видимому, быстрее в C (т. lapply использовании функции). Хотя разница кажется большой, мы все же остаемся на уровне микросекунды, а дорогостоящей является оценка функции R на каждой итерации. Простой пример:

   ffR = function(x)  {
       ans = vector("list", length(x))
       for(i in seq_along(x)) ans[[i]] = x[[i]]
       ans 
   }
   
   ffC = inline::cfunction(sig = c(R_x = "data.frame"), body = '
       SEXP ans;
       PROTECT(ans = allocVector(VECSXP, LENGTH(R_x)));
       for(int i = 0; i < LENGTH(R_x); i++) 
              SET_VECTOR_ELT(ans, i, VECTOR_ELT(R_x, i));
       UNPROTECT(1);
       return(ans); 
   ')
   
   set.seed(007) 
   myls = replicate(1e3, runif(1e3), simplify = FALSE)     
   mydf = as.data.frame(myls)
   
   all.equal(ffR(myls), ffC(myls))
   #[1] TRUE 
   all.equal(ffR(mydf), ffC(mydf))
   #[1] TRUE
   
   microbenchmark::microbenchmark(ffR(myls), ffC(myls), 
                                  ffR(mydf), ffC(mydf),
                                  times = 30)
   #Unit: microseconds
   #      expr       min        lq    median        uq       max neval
   # ffR(myls)  3933.764  3975.076  4073.540  5121.045 32956.580    30
   # ffC(myls)    12.553    12.934    16.695    18.210    19.481    30
   # ffR(mydf) 14799.340 15095.677 15661.889 16129.689 18439.908    30
   # ffC(mydf)    12.599    13.068    15.835    18.402    20.509    30
   

2. Как уже упоминалось @Roland, он запускает скомпилированный цикл C, а интерпретируемый цикл R

Хотя при векторизации кода вы должны учитывать некоторые вещи.

1. Если ваш набор данных (пусть его называют df) имеет класс data.frame, некоторые векторизованные функции (такие как colMeans, colSums, rowSums и т.д.) colSums rowSums преобразовать его в матрицу, просто потому, что именно так они были разработаны, Это означает, что для большого df это может создать огромные накладные расходы. В то время как lapply не придется делать это, поскольку он извлекает фактические векторы из df (поскольку data.frame - это всего лишь список векторов), и, таким образом, если у вас не так много столбцов, но много строк, lapply(df, mean) иногда может быть лучшим вариантом, чем colMeans(df).
2. Еще одна вещь, которую следует помнить, состоит в том, что R имеет большое количество различных типов функций, таких как .Primitive и generic (S3, S4), см. Здесь дополнительную информацию. Общая функция должна выполнять отправку метода, которая иногда является дорогостоящей операцией. Например, mean - это общая функция S3 тогда как sum является Primitive. Таким образом, несколько раз lapply(df, sum) может быть очень эффективным по сравнению с colSums из перечисленных выше причин