如何以快速高效的方式在矩阵的每一行中将最小值赋为0？

Question

有一个矩阵 Lambda，它有 p 列和 n 行，对于每一行都希望将所有的值赋为0，除了第一列的值和其他列中最大的值（在避免第一列后的所有 p - 2 个最小值）。

目前我正在使用 for 循环来完成这个操作，如下所示：

set.seed(60)
(Lambda = matrix(sample.int(30),5))
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]   19   20   27   18   15   25
[2,]   16   28    1    4   22    7
[3,]    2   10    8   23    3   12
[4,]    5    6    9   17   11   29
[5,]   26   30   24   13   14   21

m <- ncol(Lambda) - 2
for(ir in seq_len(nrow(Lambda))){
    Lambda[ir, match(tail(sort(abs(Lambda[ir, 2:ncol(Lambda)]), decreasing = TRUE), m), Lambda[ir,])] <- 0
}
Lambda
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]   19    0   27    0    0    0
[2,]   16   28    0    0    0    0
[3,]    2    0    0   23    0    0
[4,]    5    0    0    0    0   29
[5,]   26   30    0    0    0    0

好的，一行可以得到目标，但如果有很多行，就会出现问题。有没有不使用for循环的解决方案？它可能使用lapply，但我不确定这是否真的有效。也许可以在将矩阵转换后使用data.table？

谢谢！

- iago

---

你要使用的实际矩阵有多大？ - s_baldur

实际上，我的矩阵并不是很大（无论如何，少于100万行）。这对我来说不是问题，但我更希望以更好的方式解决这个问题。 - iago

Answer 1

关于这个问题

有没有不使用for循环的解决方案。

对于某些算法，你必须编写一个for循环。这没关系！将for循环换成lapply之类的东西并不能真正提高性能（请参见https://dev59.com/QlgQ5IYBdhLWcg3wqF2l#42440872）。

然而，可以加速你的代码：

# your example
set.seed(60)
Lambda = matrix(sample.int(30),5)

original <- function(Lambda) {
  m <- ncol(Lambda) - 2
  for (ir in seq_len(nrow(Lambda))){
    Lambda[ir, match(tail(sort(abs(Lambda[ir, 2:ncol(Lambda)]), decreasing = TRUE), m), Lambda[ir,])] <- 0
  }
  Lambda
}
original(Lambda)

# a faster alternative
proposal <- function(Lambda) {

  nc <- ncol(Lambda)
  for (i in seq_len(nrow(Lambda))) {
    m <- which.max(abs(Lambda[i, -1L]))
    Lambda[i, (2:nc)[-m]] <- 0
  }
  Lambda
}
proposal(Lambda)

让我们对这两种方法进行基准测试：

bch <- bench::mark(
  original(Lambda),
  proposal(Lambda)
)
summary(bch, relative = TRUE)
# A tibble: 2 x 13
  expression              min median `itr/sec` mem_alloc `gc/sec` n_itr  n_gc 
  <bch:expr>            <dbl>  <dbl>     <dbl>     <dbl>    <dbl> <int> <dbl> 
1 original(Lambda)       25.7   24.1       1           1     1     1447     4
2 proposal(Lambda)        1      1        23.6         1     2.57  9997     3

所以，proposal 比你的原始解决方案快了约24倍（原始方案的中位时间为313.8微秒，而 proposal 为13.1 微秒）。如果这还不够快，或许可以寻找用C或C++实现了此功能的软件包。我曾尝试使用matrixStats，但没有成功。另外，你还可以使用Rcpp 将其移植到C++中，这也应该会加速代码。

- Vandenman

---

非常感谢！这已经是一个很大的改进了。 - iago

Answer 2

这里有一种选项，似乎比在 600k 行上使用 proposal() 更快，但速度提升不超过 15%：

foo <- function(Lambda) {
  nr <- nrow(Lambda)
  keep <- c(seq_len(nr), apply(Lambda[, -1], 1, which.max)*nr + seq_len(nr))
  replace(Lambda, -keep, 0L)
}

编辑

一个巨大的改进是使用markus建议的max.col()替换apply()+which.max()组合：

foo2 <- function(Lambda) {
  nr <- nrow(Lambda)
  keep <- c(seq_len(nr), max.col(Lambda[, -1], ties.method = "first")*nr + seq_len(nr))
  replace(Lambda, -keep, 0L)  
}

（更新）基准测试：

set.seed(60)
Lambda = matrix(sample.int(36e5), ncol = 6)
bench::mark(
  foo(Lambda),
  proposal(Lambda),
  foo2(Lambda),
  relative = TRUE
)[1:5]

  expression         min median `itr/sec` mem_alloc
  <bch:expr>       <dbl>  <dbl>     <dbl>     <dbl>
1 foo(Lambda)       17.7   12.1      1.09      3.67
2 proposal(Lambda)  19.3   13.1      1         1   
3 foo2(Lambda)       1      1       13.5       3.75

- s_baldur

---

如果您使用max.col(Lambda[,-1], ties.method = "first")而不是apply(Lambda[, -1], 1, which.max)，您已经很好的方法可以进一步改进。 - markus

@markus 哇，差别真是天壤之别。谢谢你，我都忘了 max.col() 了。 - s_baldur

非常感谢sindri_baldur和@markus。 - iago

Answer 3

一种可以思考的向量化版本：

vectorized <- function(Lambda) {
  max <- matrixStats::rowMaxs(Lambda, cols = -1)
  noreplace <- sweep(Lambda, 1, max, "==")
  noreplace[, 1] <- TRUE
  Lambda * noreplace
}

但它并不比 @Vandenman 的 for 循环更快。

- F. Privé

---

非常感谢！它运行得很好。实际上，与sindri_baldur或更大的基准测试大小相比，您的函数几乎与他和markus的foo2一样快。 - iago