r - 当数据在列表中时，我可以向量化代码吗？

我正在优化代码，但遇到了一些问题。我知道R中最大的提速来自矢量化代码，而不是使用循环。但是，我的数据在列表中，并且不确定是否可以对代码进行矢量化处理。我尝试使用apply函数(例如lapply，vapply)，但我读到这些函数仅用于编写更简洁的代码，并且实际上是在后台使用循环!

这是我代码中的三个最大瓶颈，尽管我认为第一部分什么也做不了。

1)读取数据

我可以批量处理尺寸为277x349的1000个矩阵。这是我脚本中最大的瓶颈，但是我使用doMC包通过foreach函数利用了多个内核，从而稍微缓解了这个问题。这将导致包含1000 277x349矩阵的列表。

出于这个问题的目的，假设我们有1000个尺寸为277 x 349的矩阵的列表

# Fake data
l <- list()
for(i in 1:1000) {
  l[[i]] <- matrix(rnorm(277*349), nrow=277, ncol=349)
}

2)瓶颈1

我需要比较一些引用矩阵(相同尺寸)。这导致将我的列表中的1000个矩阵与引用矩阵进行比较，以获得1000个距离的向量。如果我知道矩阵的维数相同，是否可以对这一步骤进行矢量化处理？

这是一些代码:

# The reference matrix
r <- matrix(rnorm(277*349), nrow=277, ncol=349)
# The number of non NA values in matrix. Do not need to worry about this...
K <- 277*349

# Make a function to calculate distances
distance <- function(xi, xj, K, na.rm=TRUE) {
  sqrt(sum((xi - xj)^2, na.rm=na.rm)/K)
}

# Get a vector containing all the distances
d <- vapply(l, distance, c(0), xj=r, K=K)

使用vapply可以非常快地完成此步骤，但这是代码的第三慢的部分。

3)瓶颈2

现在，我想对我的引用矩阵进行J个“最接近”的矩阵的加权平均矩阵。 (有一个排序步骤，但为简单起见，假设使用d[1] < d[2] < ... < d[1000])。我想获得当J = 1,2，...，1000时的加权平均矩阵

# Get the weighted matrix
weightedMatrix <- function(listOfData, distances, J) {
  # Calculate weights:
  w <- d[1:J]^{-2} / sum(d[1:J]^{-2})

  # Get the weighted average matrix
  # *** I use a loop here ***
  x_bar <- matrix(0, nrow=nrow(listOfData[[1]]), ncol=ncol(listOfData[[1]]))
  for(i in 1:J) {
    x_bar <- x_bar + {listOfData[[i]] * w[i]}
  }

  return(x_bar)
}

# Oh no! Another loop...
res <- list()
for(i in 1:length(l) ) {
  res[[i]] <- weightedMatrix(l, d, J=i)
}

我有点难过。我没有看到一种对矩阵列表进行矢量化处理的直观方法。

我正在编写的脚本会经常被调用，因此即使有一点改进也可以加起来!

编辑:

RE:1)读取数据

我忘了提到我的数据采用特殊格式，因此我必须使用特殊的数据读取功能来读取R中的数据。文件采用netcdf4格式，并且我正在使用nc_open包中的ncdf4函数进行访问文件，然后我必须使用ncvar_get函数来读取感兴趣的变量。不错的是，可以从磁盘读取文件中的数据，然后我可以使用ncvar_get将数据读取到内存中，以使用R对它们进行操作。

话虽这么说，尽管我知道矩阵的大小以及将要拥有的矩阵个数，但我还是用数据列表问了我的问题，因为foreach函数使我能够进行并行计算，可以输出并行化循环的结果在列表中。我发现使用foreach函数，数据读取步骤快了大约3倍。

我想我以后可以将数据安排为3d数组，但是也许分配3d数组所花费的时间可能比节省的时间更长？我明天将不得不尝试。

编辑2:

以下是一些我编写脚本的时间。

原始脚本:

[1] "Reading data to memory"
user  system elapsed
176.063  44.070  26.611

[1] "Calculating Distances"
user  system elapsed
2.312   0.000   2.308

[1] "Calculating the best 333 weighted matrices"
user  system elapsed
63.697  28.495   9.092

到目前为止，我进行了以下改进:(1)根据Martin Morgan的建议，在读取数据之前预先分配了列表，(2)改进了加权矩阵计算。

[1] "Reading data to memory"
user  system elapsed
192.448  38.578  27.872

[1] "Calculating Distances"
user  system elapsed
2.324   0.000   2.326

[1] "Calculating all 1000 weighted matrices"
user  system elapsed
1.376   0.000   1.374

一些注意事项:

我在foreach循环中使用12个内核来读取数据(registerDoMC(12))。整个脚本在改进前后需要大约40s/36s的时间来运行。

我的瓶颈2的时机已大大改善。以前，我只计算加权矩阵的前三分之一(即333)，但是现在脚本可以只计算原始时间的一小部分中的所有加权矩阵。

感谢您的帮助，稍后我将尝试调整代码，以查看是否可以更改脚本以使用3D数组而不是列表。我现在要花一些时间来验证计算结果，以确保它们可以正常工作!

最佳答案

我的“低落果实”(scan；预分配和填充)似乎无关紧要，所以...

距离计算中的操作看起来对我来说足够矢量化了。也许您可以对所有矩阵进行一次距离计算，从而挤出一些额外的速度，但这可能会使代码难以理解。

weightedMatrix计算看起来还有改进的余地。让我们计算一下

w <- d^(-2) / cumsum(d^(-2))

对于加权矩阵m，我认为连续矩阵之间的关系只是m' = m * (1 - w[i]) + l[[i]] * w[i]，所以

res <- vector("list", length(l))
for (i in seq_along(l))
    if (i == 1L) {
        res[[i]] = l[[i]] * w[[i]]
    } else  {
        res[[i]] = res[[i - 1]] * (1 - w[[i]])  + l[[i]] * w[[i]]
    }

这将res的计算从二次更改为线性。我关于优于线性性能的想法只是(可能还误导了)预感。我还没有追求。

返回到预分配和填充以及@flodel的注释，我们有

f0 <- function(n) {
    ## good: pre-allocate and fill
    l = vector("list", n)
    for (i in seq_along(l))
        l[[i]] = 1
    l
}

f1 <- function(n) {
    ## bad: copy and append
    l = list()
    for (i in seq_len(n))
        l[[i]] = 1
    l
}

产生相同的结果

> identical(f0(100), f1(100))
[1] TRUE

但是有不同的表现

> sapply(10^(1:5), function(i) system.time(f0(i))[3])
elapsed elapsed elapsed elapsed elapsed
  0.000   0.000   0.002   0.014   0.134
> sapply(10^(1:5), function(i) system.time(f1(i))[3])
elapsed elapsed elapsed elapsed elapsed
  0.000   0.001   0.005   0.253  24.520

即使这对当前问题的规模无关紧要，也无关紧要，但似乎人们应该采用更好的预分配和填充策略，以避免不得不猜测它是否相关。更好的是，使用*apply或在这种情况下使用replicate系列以避免考虑它

l <- replicate(1000, matrix(rnorm(277*349), nrow=277, ncol=349), simplify=FALSE)

关于r - 当数据在列表中时，我可以向量化代码吗？，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/17077183/