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

Question

也许这个问题非常愚蠢.

我试图"矢量化"以下循环:

set.seed(0)
x <- 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)) x[i] <- x[sig[i]]
x
# [1] 0.90 0.27 0.66 0.91 0.66 0.91 0.94 0.91 0.94 0.63

我认为这很简单,x[sig]但结果并不匹配.

set.seed(0)
x <- round(runif(10), 2)
x[] <- x[sig]
x
# [1] 0.90 0.27 0.66 0.91 0.37 0.57 0.94 0.20 0.90 0.63

怎么了？

---

备注

显然,从输出中我们看到for循环并且x[sig]是不同的.后者的含义很明确:排列,因此许多人倾向于认为循环只是做了一些错误的东西.但永远不要那么肯定; 它可以是一些定义明确的动态过程.本问答的目的不是判断哪个是正确的,而是解释为什么它们不相同.希望它为理解"矢量化"提供了坚实的案例研究.

Answer 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]在第一条指令中进行修改,然后将该修改后的值传递给x[3]第二条指令.所以x[3]会有相同的价值x[1].但是,在并行执行中,x[3]等于x[2].结果,这个循环不能"矢量化".

在"矢量化"理论中,

• 示例1在数据中具有"写后读"依赖性: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,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]

在第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>"

我R对话已经被分配地址所指向的存储块<0x28f7340>的x,你可能会看到一个不同的值,当你运行代码.但是,tracemem循环后输出不会改变,这意味着没有副本x.所以循环确实在不使用额外内存的情况下进行"就地"修改.

但是,循环没有进行"就地"排列."就地"排列是一种更复杂的操作.不仅x需要在循环中交换元素,sig还需要交换元素(最后sig也是如此1:10).