在R中高效地进行矩阵运算

Question

以下代码是我项目的核心，但是考虑到我的问题规模，目前它太慢了。是否有更有效的方法来实现相同的结果？

nbassets <- 80
nbrisksource <- 100
nbsimul <- 300000
set.seed(100)
#generate random number for each 100 source of risk in many simulations
random <- matrix(runif(nbsimul*nbrisksource)+0.9,nrow=nbsimul,ncol=nbrisksource)
# random vulnerability to each source of risk for each of 120 assets
EL_decomp <- matrix(runif(nbassets*nbrisksource),nrow=nbassets,ncol=nbrisksource)
#initiate matrix to store asset returns
asset_ret <- matrix(NA, nrow=nbsimul,ncol=nbassets)


ptm <- proc.time()
#loop through each asset
 for (i in 1:nbassets){
  #determine if the asset has been impacted by any source of risk, if yes return is -1, otherwise 0
  asset_ret[,i] <- apply(matrix(EL_decomp[i,], nrow=nbsimul,ncol=nbrisksource,byrow=TRUE) < random,1,all)-1
}
print(proc.time() - ptm)

ptm <- proc.time()

- Simon Wuya

Answer 1

我可以将其加速18倍，基本上跳过所有矩阵写入并利用R的循环机制：

n_80     <- 80
n_100    <- 100
n_300000 <- 300000
set.seed(100)
mat_300000_100 <- matrix(runif(n_300000*n_100), nrow=n_300000, ncol=n_100)
mat_80_100     <- matrix(runif(n_80    *n_100), nrow=n_80,     ncol=n_100)
mat_300000_80  <- matrix(NA, nrow=n_300000, ncol=n_80)

首先，要删除矩阵，因为>将向量进行循环利用。必须转置，因为>按列而不是按行应用向量。如果可能的话，请使用优化的函数，例如colSums而不是apply。这里可以将apply(v,2,'all')替换为colSums(v)==length_v。

ptm <- proc.time()
for (i in 1:n_80) mat_300000_80[,i] <- colSums(mat_80_100[i,] < t(mat_300000_100))==n_100-1
print(proc.time() - ptm) # 17s

最终，在循环外只需执行一次transpose()。（或者由于你的示例中的值是完全随机的，甚至可以不执行...）

ptm <- proc.time()
mat_100_300000 <- t(mat_300000_100)
for (i in 1:n_80) mat_300000_80[,i] <- colSums(mat_80_100[i,] < mat_100_300000)==n_100-1
print(proc.time() - ptm) # 8s

- Arthur

---

你是对的！但即使使用赋值，apply 的性能提升也不到0.5秒。 - Arthur

谢谢。尽管您的代码速度更快，但它并没有执行相同的操作。也就是说，在我的代码中，all()函数是逐行应用的，而您则将其应用于整个矩阵。结果并不相同。虽然我选择的示例不太好，但我已经进行了编辑，以便清楚地表明这两个操作会导致不同的答案。 - Simon Wuya

只是为了明确起见：apply(matrix(EL_decomp[i,],nrow=nbsimul,ncol=nbrisksource,byrow=TRUE) 返回一个布尔数组，而 all(mat_80_100[i,] < mat_100_300000) 只返回一个布尔值。 - Simon Wuya

@TheTime，你差点就做到了，但是在你优雅的一行代码中有一个错误（将nbassets替换为nbrisksource），像这样：apply(EL_decomp, 1, function(x) (colSums(x < t(random)) == nbrisksource) - 1L)。如果你不介意把这个作为答案发出来，我会接受它。否则我会自己发布。再次感谢。 - Simon Wuya

@SimonWuya，你说得对，我没有看到那个微妙之处。已经更正了，而且速度还很快。 - Arthur

Answer 2

事情可以得到极大的改善。以下是旧代码和新代码的比较：

nbassets     <- 80
nbrisksource    <- 100
nbsimul <- 300000
set.seed(100)
random <- matrix(runif(nbsimul*nbrisksource)+0.9, nrow=nbsimul,ncol=nbrisksource)

EL_decomp     <- matrix(runif(nbassets    *nbrisksource), nrow=nbassets,     ncol=nbrisksource)
asset_ret1  <- matrix(NA, nrow=nbsimul, ncol=nbassets)
asset_ret2  <- matrix(NA, nrow=nbsimul, ncol=nbassets)

ptm <- proc.time()
for (i in 1:nbassets){
  #determine if the asset has been impacted by any source of risk, if yes return is -1, otherwise 0
  asset_ret1[,i] <- apply(matrix(EL_decomp[i,],nrow=nbsimul,ncol=nbrisksource,byrow=TRUE) < random,1,all)-1
}
print(head(asset_ret1))
print(proc.time() - ptm) #182s on my old mac

#improved version
ptm <- proc.time()
randomt <- t(random)
asset_ret2 <- apply(EL_decomp, 1, function(x) (colSums(x < randomt) == nbrisksource))- 1L
print(head(asset_ret2))
print(proc.time() - ptm) #14s
print(identical(asset_ret1,asset_ret2)) 

- Simon Wuya

---

请注意：在出现NA的情况下，apply(x, 1, all)和colSums(x) == n之间可能会存在不匹配的情况；您可能需要先计算每列的非NA值，然后再假设它是“n”。 - alexis_laz