r - 为什么“矢量化”这个简单的 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

“矢量化”在第一个示例中成功，但在第二个示例中不成功。为什么？

这里是谨慎的分析。“向量化”从循环展开开始，然后并行执行几条指令。循环是否可以“向量化”取决于循环所携带的数据依赖关系。

展开示例 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]第 2 条指令。所以x[3]将具有与相同的值x[1]。但是，在并行执行中，x[3]等于x[2]. 结果，这个循环不能被“矢量化”。

在“矢量化”理论中，

示例 1 在 data 中有一个“write-after-read”依赖：在x[i]读取后修改；
示例 2 在数据中具有“写后读”依赖性：x[i]在修改后读取。

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

深入

可能很多人到现在都糊涂了。“向量化”是“并行处理”吗？

是的。在 1960 年代，当人们想知道应该为高性能计算设计什么样的并行处理计算机时，Flynn 将设计思想分为 4 类。“SIMD”（单指令，多数据）类别被称为“矢量化”，具有“SIMD”功能的计算机称为“矢量处理器”或“阵列处理器”。

在 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。“向量化”不再明确指代“SIMD”。R不是一种我们可以对 CPU 寄存器进行编程的语言。R中的“矢量化”只是对“SIMD”的类比。在之前的问答中：“向量化”一词在不同的上下文中是否意味着不同的东西？我试图解释这一点。下面的地图说明了这个类比是如何产生的：

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

除了“数据依赖”，还有“控制依赖”，即在“向量化”中处理“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]

第 3 条和第 5 条指令之间存在“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 中实现循环被视为“矢量化”，那就顺其自然吧。但是没有机会用“SIMD”进一步“矢量化”C/C++循环。

备注 2

这个问答是由这个问答推动的。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列，而要查找的列在第一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>我的 R 会话分配了一个由地址指向的内存块，x运行代码时您可能会看到不同的值。但是，tracemem循环后的输出不会改变，这意味着没有复制x。所以循环确实在不使用额外内存的情况下进行“就地”修改。

但是，循环并没有进行“就地”排列。“就地”置换是一个更复杂的操作。不仅x需要沿着循环交换元素，sig还需要交换元素（最后sig是1:10）。

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

问题描述

解决方案

暖身

深入

一张完整的图片

回到问题中的例子

备注 1

备注 2

备注 3

推荐阅读