为什么“矢量化”这个简单的R循环会产生不同的结果?

热身

作为热身,考虑两个简单的例子。

## 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

这是审慎的分析。”“矢量化”从循环展开开始,然后并行执行多条指令。循环是否可以“矢量化”取决于循环所承载的数据依赖关系。

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] 等于 . 因此,此循环无法“矢量化”。

在“矢量化”理论中,

• 示例1在数据中具有“先读后写”的依赖关系: x[i] 读取后修改;
• 示例2在数据中具有“先读后写”的依赖关系: x[i] 修改后读取。

具有“先读后写”数据依赖关系的循环可以“矢量化”,而具有“先读后写”数据依赖关系的循环则不能。

深入

也许现在很多人都困惑了。”“矢量化”是“并行处理”吗?

在1960年代,编程语言并不多。人们编写汇编程序(当编译器发明时是FORTRAN)直接对CPU寄存器进行编程。“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 R 不是一种可以编程CPU寄存器的语言。R中的“矢量化”只是“SIMD”的一个类比。在以前的问答中;答: Does the term "vectorization" mean different things in different contexts?

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

所以,示例1中循环的R“矢量化”类似于

## 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 在内存中没有混叠,即存储向量的内存块 b 和 不要重叠。这是一个理想的情况,因为没有内存别名意味着没有数据依赖性。

回到问题中的例子

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]

第3条和第5条指令之间存在“先读后写”的数据依赖关系,因此循环无法“矢量化”(请参阅) 备注1

那么,那是什么呢 x[] <- x[sig]

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

自 "["

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[信号] 相当于

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

备注1

this Q & A . OP最初提出了一个循环

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"]]

但这可能是错误的。后来OP把循环改成

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

这消除了内存混叠问题,因为要替换的列是第一列 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>"

<0x28f7340> 对于 十 运行代码时,可能会看到不同的值。但是 循环后不会改变,这意味着 是制造的。因此,循环确实在不使用额外内存的情况下进行“就地”修改。

但是,循环没有进行“就地”排列“就地”置换是一个更复杂的操作。不仅仅是 十 sig 也需要交换(最后, 信号 会是 1:10 ).

有一个更简单的解释。使用循环,您将覆盖 x sample(x, replace=TRUE) )--你是否需要这样一个复杂的问题,取决于你想要实现什么。

对于矢量化代码,您只需要对 十 (没有替换),这就是你得到的。矢量化代码是 不 十 :

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

这里没有混淆的危险: x2 最初引用同一个内存位置,但一旦您更改 x2个

这与内存块别名无关(我以前从未遇到过这个术语)。以一个特定的排列为例,遍历在C语言或汇编语言(或任何语言)级别上的实现都会发生的赋值;它是任何顺序for循环的行为与任何“真正”排列的行为之间的内在联系 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 

有趣的是,尽管R“矢量化”不同于“SIMD”(正如OP很好地解释的),但在确定循环是否“可矢量化”时,同样的逻辑也适用。下面是一个演示,使用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

不 “可矢量化”。

// "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