r - 使用 for 循环进行矩阵计算

Question

我是 R 的新手，由于许可证验证问题而从 GAUSS 迁移。

我想加速以下创建n×k矩阵A的代码。给定n×1向量x和参数向量mu、sig(均为k维)， A 创建为 A[i,j]=dnorm(x[i], mu[j], sigma[j])。以下代码对于较小的数字 n=40, k=4 可以正常工作，但当 n 约为 10^6 和 时，速度会显着减慢>k 的大小与 n^{1/3} 大致相同。

我正在做模拟实验来验证bootstrap的有效性，所以我需要重复计算矩阵A#ofsimulation × #bootstrap次，并且花费的时间很少因为我想尝试许多不同的 n,k 值。我尽可能多地对代码进行矢量化(感谢 dnorm 的矢量参数)，但是我可以要求更快吗？

先行感谢您的帮助。

x   = rnorm(40)
mu  = c(-1,0,4,5)
sig = c(2^2,0.5^2,2^2,3^2)
n   = length(x)
k   = length(mu)    
A   = matrix(NA,n,k)

for(j in 1:k){
    A[,j]=dnorm(x,mu[j],sig[j])
}

Answer 1

您的方法可以放入这样的函数中

A.fill <- function(x,mu,sig) {
  k <- length(mu)   
  n <- length(x)
  A <- matrix(NA,n,k)
  for(j in 1:k) A[,j] <- dnorm(x,mu[j],sig[j])
  A

}

很明显，您正在填充矩阵 A逐列。 R按列存储矩阵的条目(就像 Fortran)。 这意味着矩阵可以通过一次调用 dnorm 来填充。使用 x 的适当重复, mu ，和sig 。矢量z将堆叠所需矩阵的列。然后只需指定行数和列数即可从该向量形成要返回的矩阵。看下面的函数

B.fill <- function(x,mu,sig) { 
  k <- length(mu)
  n <- length(x)
  z <- dnorm(rep(x,times=k),rep(mu,each=n),rep(sig,each=n))
  B <- matrix(z,nrow=n,ncol=k)
  B

}

让我们用您的数据制作一个示例并进行测试，如下所示:

N <- 40
set.seed(11)
x <- rnorm(N)
mu <- c(-1,0,4,5)
sig <- c(2^2,0.5^2,2^2,3^2)
A <- A.fill(x,mu,sig)
B <- B.fill(x,mu,sig)

all.equal(A,B)

# [1] TRUE

我假设n是 k 的整数倍.

添加

正如评论中所指出的，对于 n 的大值，B.fill 相当慢。 。 原因在于构造 rep(...,each=...) .

那么有没有办法提速A.fill 。 我测试了这个功能:

C.fill <- function(x,mu,sig) {
  k <- length(mu)
  n <- length(x)
  sapply(1:k,function(j) dnorm(x,mu[j],sig[j]), simplify=TRUE)
}

该函数比 A.fill 快约 20%。