【问题标题】：Vectorize min() for matrix为矩阵向量化 min()
【发布时间】：2016-05-19 18:28:50
【问题描述】：

我希望矢量化以下循环：

for (i in 1:n) {
    for (j in 1:m) {
        temp_mat[i,j]=min(temp_mat[i,j],1);
    }
}

我以为我可以做到temp_mat=min(temp_mat,1)，但这并没有给我想要的结果。有没有办法对这个循环进行矢量化以使其更快？

【问题讨论】：

标签： r matrix vectorization

【解决方案1】：

只需使用temp_mat <- pmin(temp_mat, 1)。有关并行最小值的更多使用，请参阅?pmin。

例子：

set.seed(0); A <- matrix(sample(1:3, 25, replace = T), 5)
#> A
#     [,1] [,2] [,3] [,4] [,5]
#[1,]    3    1    1    3    3
#[2,]    1    3    1    2    3
#[3,]    2    3    1    3    1
#[4,]    2    2    3    3    2
#[5,]    3    2    2    2    1
B <- pmin(A, 2)
#> B
#     [,1] [,2] [,3] [,4] [,5]
#[1,]    2    1    1    2    2
#[2,]    1    2    1    2    2
#[3,]    2    2    1    2    1
#[4,]    2    2    2    2    2
#[5,]    2    2    2    2    1

更新

由于您有计算科学背景，我想提供更多信息。

pmin 速度很快，但远非高性能。它的前缀“parallel”仅暗示element-wise。 R中“向量化”的含义与HPC中的“SIMD向量化”不同。 R 是一种解释型语言，所以 R 中的“向量化”意味着选择 C 级循环而不是 R 级循环。因此，pmin 只是用一个普通的 C 循环编码。

真正的高性能计算应该受益于 SIMD 矢量化。我相信你知道 SSE/AVX 内在函数。因此，如果您编写一个简单的 C 代码，使用来自SSE2 的_mm_min_pd，您将获得从pmin 提高约2 倍的速度；如果您从 AVX 看到 _mm256_min_pd，您将从 pmin 获得约 4 倍的加速。

很遗憾，R 本身不能执行任何 SIMD。 我在Does R leverage SIMD when doing vectorized calculations? 上就这个问题回复了帖子。对于您的问题，即使您将 R 链接到 HPC BLAS，pmin 也不会从 SIMD 中受益，因为pmin 不涉及任何 BLAS 操作。所以更好的选择是自己编写编译后的代码。

【讨论】：

我认为这与 R 是主要列无关。 pmin 返回传递的项目之间的元素最小值，因此 1 只是被强制转换为相同尺寸的对象。例如，如果传递了不同大小的向量，则会收到有关元素被回收的警告。
虽然这不是 OP 所要求的，但有许多方法可以在 R 中实现大规模并行计算。示例包括包snow、multicore、foreach、Rmpi 、Rth、gputools 和基本包parallel，它合并了其中一些以前独立的包。这是一个good reference，关于该主题的最新技术。
"例如，OP 的问题不能从这种形式中受益，变成它受内存限制。设置多个线程会导致速度变慢。"为什么？将矩阵拆分为单独的块并让每个块由单独的线程处理不是并行编程中的标准任务吗？在这种情况下，线程之间不需要通信。在我看来，这就像一个令人尴尬的并行情况的标准示例，可以用 R 轻松处理。
没问题。如果你愿意，我们可以把它留在那里。我不想引起争议，我只是不明白你在这种情况下的观点。与您提到的 Cholesky 分解相比，这在我看来是完全不同的情况。

【解决方案2】：

这个问题有点令人困惑，因为min() 是矢量化的。但是，为了在这种特定情况下获得所需的结果，不必使用此功能。 逻辑子集提供了一种可能更有效（当然也更紧凑）的替代方案。

如果我正确理解了您想要的输出，则可以使用单个命令来修改代码中嵌套循环所执行的矩阵：

temp_mat[temp_mat > 1] <- 1

希望这会有所帮助。

set.seed(123)
temp_mat <- matrix(2*runif(50),10)
#> temp_mat
#           [,1]       [,2]      [,3]       [,4]      [,5]
# [1,] 0.5751550 1.91366669 1.7790786 1.92604847 0.2856000
# [2,] 1.5766103 0.90666831 1.3856068 1.80459809 0.8290927
# [3,] 0.8179538 1.35514127 1.2810136 1.38141056 0.8274487
# [4,] 1.7660348 1.14526680 1.9885396 1.59093484 0.7376909
# [5,] 1.8809346 0.20584937 1.3114116 0.04922737 0.3048895
# [6,] 0.0911130 1.79964994 1.4170609 0.95559194 0.2776121
# [7,] 1.0562110 0.49217547 1.0881320 1.51691908 0.4660682
# [8,] 1.7848381 0.08411907 1.1882840 0.43281587 0.9319249
# [9,] 1.1028700 0.65584144 0.5783195 0.63636202 0.5319453
#[10,] 0.9132295 1.90900730 0.2942273 0.46325157 1.7156554
temp_mat[temp_mat > 1] <- 1
#> temp_mat
#           [,1]       [,2]      [,3]       [,4]      [,5]
# [1,] 0.5751550 1.00000000 1.0000000 1.00000000 0.2856000
# [2,] 1.0000000 0.90666831 1.0000000 1.00000000 0.8290927
# [3,] 0.8179538 1.00000000 1.0000000 1.00000000 0.8274487
# [4,] 1.0000000 1.00000000 1.0000000 1.00000000 0.7376909
# [5,] 1.0000000 0.20584937 1.0000000 0.04922737 0.3048895
# [6,] 0.0911130 1.00000000 1.0000000 0.95559194 0.2776121
# [7,] 1.0000000 0.49217547 1.0000000 1.00000000 0.4660682
# [8,] 1.0000000 0.08411907 1.0000000 0.43281587 0.9319249
# [9,] 1.0000000 0.65584144 0.5783195 0.63636202 0.5319453
#[10,] 0.9132295 1.00000000 0.2942273 0.46325157 1.0000000

【讨论】：