Почему "векторизация" этой простой циклы R дает другой результат?

Ответ 1

разогреть

В качестве разминки рассмотрим два более простых примера.

## example 1
x <- 1:11
for (i in 1:10) x[i] <- x[i + 1]
x
# [1]  2  3  4  5  6  7  8  9 10 11 11

x <- 1:11
x[1:10] <- x[2:11]
x
# [1]  2  3  4  5  6  7  8  9 10 11 11

## example 2
x <- 1:11
for (i in 1:10) x[i + 1] <- x[i]
x
# [1] 1 1 1 1 1 1 1 1 1 1 1

x <- 1:11
x[2:11] <- x[1:10]
x
# [1]  1  1  2  3  4  5  6  7  8  9 10

"Векторизация" успешна в первом примере, но не во втором. Зачем?

Вот осмотрительный анализ. "Векторизация" начинается с разворачивания цикла, а затем выполняет несколько инструкций параллельно. Независимо от того, может ли цикл быть "векторизован", зависит от зависимостей данных, переносимых циклом.

Развертывание цикла в примере 1 дает

x[1]  <- x[2]
x[2]  <- x[3]
x[3]  <- x[4]
x[4]  <- x[5]
x[5]  <- x[6]
x[6]  <- x[7]
x[7]  <- x[8]
x[8]  <- x[9]
x[9]  <- x[10]
x[10] <- x[11]

Выполнение этих инструкций один за другим и их выполнение одновременно дают идентичный результат. Таким образом, этот цикл может быть "векторизован".

Цикл в примере 2

x[2]  <- x[1]
x[3]  <- x[2]
x[4]  <- x[3]
x[5]  <- x[4]
x[6]  <- x[5]
x[7]  <- x[6]
x[8]  <- x[7]
x[9]  <- x[8]
x[10] <- x[9]
x[11] <- x[10]

К сожалению, выполнение этих инструкций один за другим и их выполнение одновременно не давали бы одинакового результата. Например, при выполнении их один за другим x[2] изменяется в 1-й инструкции, тогда это модифицированное значение передается x[3] во второй инструкции. Таким образом, x[3] будет иметь такое же значение, как x[1]. Однако при параллельном выполнении x[3] равно x[2]. В результате этот цикл не может быть "векторизован".

В теории "векторизации",

• Пример 1 имеет зависимость от записи после чтения: x[i] изменяется после его чтения;
• Пример 2 имеет зависимость "после записи" в данных: x[i] считывается после его изменения.

Цикл с зависимостью "после записи" может быть "векторизован", а цикл с зависимостью "после чтения" не может быть.

глубоко

Наверное, многие люди были смущены. "Векторизация" - это "параллельная обработка"?

Да. В 1960 году, когда люди задавались вопросом, какой компьютер с параллельной обработкой должен быть разработан для высокопроизводительных вычислений, Флинн классифицировал идеи дизайна на 4 типа. Категория "SIMD" (одна инструкция, несколько данных) - это "векторизация", а компьютер с "SIMD" - это "векторный процессор" или "процессор массива".

В 1960 году было не так много языков программирования. Люди написали сборку (тогда FORTRAN, когда был изобретен компилятор), чтобы напрямую запрограммировать регистры процессора. Компьютер "SIMD" способен загружать несколько данных в векторный регистр с помощью одной команды и выполнять одну и ту же арифметику по этим данным одновременно. Таким образом, обработка данных действительно параллельна. Рассмотрим наш пример 1 снова. Предположим, что векторный регистр может содержать два векторных элемента, тогда цикл может быть выполнен с 5 итерациями, используя векторную обработку, а не 10 итераций, как при скалярной обработке.

reg <- x[2:3]  ## load vector register
x[1:2] <- reg  ## store vector register
-------------
reg <- x[4:5]  ## load vector register
x[3:4] <- reg  ## store vector register
-------------
reg <- x[6:7]  ## load vector register
x[5:6] <- reg  ## store vector register
-------------
reg <- x[8:9]  ## load vector register
x[7:8] <- reg  ## store vector register
-------------
reg <- x[10:11] ## load vector register
x[9:10] <- reg  ## store vector register

Сегодня существует много языков программирования, таких как R. "Векторизация" уже недвусмысленно ссылается на "SIMD". R не является языком, на котором мы можем программировать регистры процессора. "Векторизация" в R является просто аналогом "SIMD". В предыдущем Q & A: означает ли термин "векторизация" разные вещи в разных контекстах? Я попытался это объяснить. Следующая карта иллюстрирует, как делается эта аналогия:

single (assembly) instruction    -> single R instruction
CPU vector registers             -> temporary vectors
parallel processing in registers -> C/C++/FORTRAN loops with temporary vectors

Таким образом, R "векторизация" цикла в примере 1 - это что-то вроде

## the C-level loop is implemented by function "["
tmp <- x[2:11]  ## load data into a temporary vector
x[1:10] <- tmp  ## fill temporary vector into x

Большую часть времени мы просто делаем

x[1:10] <- x[2:10]

без явного назначения временного вектора переменной. Созданный временной блок памяти не имеет указателя на какую-либо переменную R и, следовательно, подлежит сборке мусора.

полная картина

В приведенном выше примере "векторизация" не вводится с простейшим примером. Очень часто "векторизация" вводится с чем-то вроде

a[1] <- b[1] + c[1]
a[2] <- b[2] + c[2]
a[3] <- b[3] + c[3]
a[4] <- b[4] + c[4]

где a, b и c не сглажены в памяти, то есть блоки памяти, хранящие векторы a, b и c, не перекрываются. Это идеальный случай, так как отсутствие псевдонимов памяти не подразумевает зависимости данных.

Помимо "зависимости данных", существует также "управляющая зависимость", то есть дело с "if... else..." в "векторизации". Однако по причине времени и пространства я не буду подробно останавливаться на этом вопросе.

вернуться к примеру в вопросе

Теперь пришло время исследовать цикл в вопросе.

set.seed(0)
x <- round(runif(10), 2)
sig <- sample.int(10)
# [1]  1  2  9  5  3  4  8  6  7 10
for (i in seq_along(sig)) x[i] <- x[sig[i]]

Развертывание цикла дает

x[1]  <- x[1]
x[2]  <- x[2]
x[3]  <- x[9]   ## 3rd instruction
x[4]  <- x[5]
x[5]  <- x[3]   ## 5th instruction
x[6]  <- x[4]
x[7]  <- x[8]
x[8]  <- x[6]
x[9]  <- x[7]
x[10] <- x[10]

Существует зависимость данных "read-after-write" между третьей и пятой инструкциями, поэтому цикл не может быть "векторизован" (см. Замечание 1).

Итак, что делает x[] <- x[sig]? Пусть сначала явно выписывает временный вектор:

tmp <- x[sig]
x[] <- tmp

Поскольку "[" вызывается дважды, на самом деле два цикла C-уровня за этим "векторизованным" кодом:

tmp[1]  <- x[1]
tmp[2]  <- x[2]
tmp[3]  <- x[9]
tmp[4]  <- x[5]
tmp[5]  <- x[3]
tmp[6]  <- x[4]
tmp[7]  <- x[8]
tmp[8]  <- x[6]
tmp[9]  <- x[7]
tmp[10] <- x[10]

x[1]  <- tmp[1]
x[2]  <- tmp[2]
x[3]  <- tmp[3]
x[4]  <- tmp[4]
x[5]  <- tmp[5]
x[6]  <- tmp[6]
x[7]  <- tmp[7]
x[8]  <- tmp[8]
x[9]  <- tmp[9]
x[10] <- tmp[10]

Итак, x[] <- x[sig] эквивалентно

for (i in 1:10) tmp[i] <- x[sig[i]]
for (i in 1:10) x[i] <- tmp[i]
rm(tmp); gc()

который вовсе не является исходным циклом, данным в вопросе.

Примечание 1

Если реализация цикла в Rcpp рассматривается как "векторизация", то пусть это будет. Но нет возможности продолжить "векторизовать" цикл C/C++ с помощью "SIMD".

Замечание 2

Этот Q & A мотивирован этим Q & A. ОП изначально представлял цикл

for (i in 1:num) {
  for (j in 1:num) {
    mat[i, j] <- mat[i, mat[j, "rm"]]
  }
}

Заманчиво "векторизовать" его как

mat[1:num, 1:num] <- mat[1:num, mat[1:num, "rm"]]

но это потенциально неправильно. Позже ОП изменил цикл на

for (i in 1:num) {
  for (j in 1:num) {
    mat[i, j] <- mat[i, 1 + num + mat[j, "rm"]]
  }
}

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

Замечание 3

Я получил несколько комментариев относительно того, является ли цикл в вопросе "модификацией x месте". Да, это. Мы можем использовать tracemem:

set.seed(0)
x <- round(runif(10), 2)
sig <- sample.int(10)
tracemem(x)
#[1] "<0x28f7340>"
for (i in seq_along(sig)) x[i] <- x[sig[i]]
tracemem(x)
#[1] "<0x28f7340>"

Мой сеанс R выделил блок памяти, указанный адресом <0x28f7340> для x и при запуске кода вы можете увидеть другое значение. Тем не менее, выход tracemem не изменится после цикла, а это означает, что не производится копия x. Таким образом, цикл действительно выполняет "на месте" модификацию без использования дополнительной памяти.

Однако цикл не выполняет перестановку "на месте". Перестановка на месте - более сложная операция. Не только элементы x нужно поменять местами по петле, элементы sig также нужно поменять местами (и, в конце концов, sig будет 1:10).

Ответ 2

Существует более простое объяснение. С вашим циклом вы переписываете один элемент x на каждом шаге, заменяя его прежнее значение одним из других элементов x. Итак, вы получаете то, о чем просили. По сути, это сложная форма выборки с заменой (sample(x, replace=TRUE)) - нужно ли вам такое усложнение, зависит от того, чего вы хотите достичь.

С вашим векторизованным кодом вы просто запрашиваете определенную перестановку x (без замены), и именно это вы получаете. Векторный код не делает то же самое, что и ваш цикл. Если вы хотите достичь того же результата с помощью цикла, вам сначала нужно сделать копию x:

set.seed(0)
x <- x2 <- round(runif(10), 2)
# [1] 0.90 0.27 0.37 0.57 0.91 0.20 0.90 0.94 0.66 0.63
sig <- sample.int(10)
# [1]  1  2  9  5  3  4  8  6  7 10
for (i in seq_along(sig)) x2[i] <- x[sig[i]]
identical(x2, x[sig])
#TRUE

Здесь нет опасности сглаживания: x и x2 первоначально ссылаются на одно и то же место памяти, но его изменение изменится, как только вы измените первый элемент x2.

Ответ 3

Это не имеет ничего общего с псевдонимом блока памяти (термин, с которым я никогда не сталкивался раньше). Возьмите конкретный пример перестановки и пройдите через назначения, которые произойдут независимо от реализации на уровне языка C или сборки (или любого другого); Это неотъемлемо от того, как будет вести себя любой последовательный цикл for-loop в сравнении с тем, как произойдет любая "истинная" перестановка (что получается с x[sig]):

sample(10)
 [1]  3  7  1  5  6  9 10  8  4  2

value at 1 goes to 3, and now there are two of those values
value at 2 goes to 7, and now there are two of those values
value at 3 (which was at 1) now goes back to 1 but the values remain unchanged

... может продолжаться, но это иллюстрирует, как это обычно не будет "истинной" перестановкой и очень необычно приведет к полному перераспределению значений. Я предполагаю, что только полностью упорядоченная перестановка (из которой, я думаю, есть только одна, т.е. 10:1), может привести к созданию нового набора x, которые были уникальными.

replicate( 100, {x <- round(runif(10), 2); 
                  sig <- sample.int(10); 
                  for (i in seq_along(sig)){ x[i] <- x[sig[i]]}; 
                  sum(duplicated(x)) } )
 #[1] 4 4 4 5 5 5 4 5 6 5 5 5 4 5 5 6 3 4 2 5 4 4 4 4 3 5 3 5 4 5 5 5 5 5 5 5 4 5 5 5 5 4
 #[43] 5 3 4 6 6 6 3 4 5 3 5 4 6 4 5 5 6 4 4 4 5 3 4 3 4 4 3 6 4 7 6 5 6 6 5 4 7 5 6 3 6 4
 #[85] 8 4 5 5 4 5 5 5 4 5 5 4 4 5 4 5

Я начал задаваться вопросом, что распределение дублирования может быть в большой серии. Выглядит довольно симметрично:

table( replicate( 1000000, {x <- round(runif(10), 5); 
                            sig <- sample.int(10); 
               for (i in seq_along(sig)){ x[i] <- x[sig[i]]}; 
                            sum(duplicated(x)) } ) )

     0      1      2      3      4      5      6      7      8 
     1    269  13113 126104 360416 360827 125707  13269    294 

Ответ 4

Интересно видеть, что, хотя "векторизация" R "отличается от" SIMD "(как хорошо объяснено ОП), та же логика применяется при определении того, является ли цикл" векционируемым ". Вот демонстрационная версия, использующая примеры в OP self-answer (с небольшой модификацией).

Пример 1 с зависимостью "после записи" является "векционируемым".

// "ex1.c"
#include <stdlib.h>
void ex1 (size_t n, size_t *x) {
  for (size_t i = 1; i < n; i++) x[i - 1] = x[i] + 1;
}

gcc -O2 -c -ftree-vectorize -fopt-info-vec ex1.c
#ex1.c:3:3: note: loop vectorized

Пример 2 с зависимостью "read-after-write" не является "векционируемым".

// "ex2.c"
#include <stdlib.h>
void ex2 (size_t n, size_t *x) {
  for (size_t i = 1; i < n; i++) x[i] = x[i - 1] + 1;
}

gcc -O2 -c -ftree-vectorize -fopt-info-vec-missed ex2.c
#ex2.c:3:3: note: not vectorized, possible dependence between data-refs
#ex2.c:3:3: note: bad data dependence

Использование C99 restrict ключевое слово намекнуть составитель не блок памяти ступенчатости между тремя массивами.

// "ex3.c"
#include <stdlib.h>
void ex3 (size_t n, size_t * restrict a, size_t * restrict b, size_t * restrict c) {
  for (size_t i = 0; i < n; i++) a[i] = b[i] + c[i];
}

gcc -O2 -c -ftree-vectorize -fopt-info-vec ex3.c
#ex3.c:3:3: note: loop vectorized