r - 重复执行矩阵乘法的有效方法

Question

我正在尝试进行矩阵乘法 S_g对于每个 i，以及每个 g 与 i。这是我到目前为止所尝试过的，但需要花费大量时间才能完成。是否有一种计算效率更高的方法可以完成完全相同的事情？

此公式中需要注意的主要内容是 S_g在矩阵乘法设置中使用 X_gamma 和 Y[,i]。 X_gamma 取决于值 g 。因此，对于每个 i，我必须执行 g矩阵乘法。

逻辑如下:

• 对于每个i，需要对每个g进行计算。然后，对于每个 g，选择 X_gamma 作为 X 的子集。以下是如何确定 X_gamma。假设 g = 3。当我们查看“set[3,]”时，我们发现 B 列是唯一一个值 != 0 的列。因此，我在 X 中选择 B 列，这将是 X_gamma。

我的主要问题是，在现实中， g = 13,000 和 i = 700 。

 library(foreach)
 library(doParallel) ## parallel backend for the foreach function
 registerDoParallel()

 T = 3
 c = 100

 X <- zoo(data.frame(A = c(0.1, 0.2, 0.3), B = c(0.4, 0.5, 0.6), C = c(0.7,0.8,0.9)),
     order.by = seq(from = as.Date("2013-01-01"), length.out = 3, by = "month")) 

 Y <- zoo(data.frame(Stock1 = rnorm(3,0,0.5), Stock2 = rnorm(3,0,0.5), Stock3 = rnorm(3,0,0.5)), 
    order.by = seq(from = as.Date("2013-01-01"), length.out = 3, by = "month"))

 l <- rep(list(0:1),ncol(X))
 set = do.call(expand.grid, l)
 colnames(set) <- colnames(X)

 I = diag(T)


 denom <- foreach(i=1:ncol(Y)) %dopar% {    
    library(zoo)
    library(stats)
    library(Matrix)
    library(base)

    result = c()
    for(g in 1:nrow(set)) {
        X_gamma = X[,which(colnames(X) %in% colnames(set[which(set[g,] != 0)]))]
        S_g = Y[,i] %*% (I - (c/(1+c))*(X_gamma %*% solve(crossprod(X_gamma)) %*% t(X_gamma))) %*% Y[,i] 
        result[g] = ((1+c)^(-sum(set[g,])/2)) * ((S_g)^(-T/2))
    }
    sum(result) 
 }

感谢您的帮助!

Answer 1

最明显的问题是您犯了一个经典错误:没有预先分配输出向量结果。对于大型向量来说，一次附加一个值可能非常低效。

在您的情况下，结果不需要是向量:您可以将结果累积为单个值:

result = 0
for(g in 1:nrow(set)) {
    # snip
    result = result + ((1+c)^(-sum(set[g,])/2)) * ((S_g)^(-T/2))
}
result

但我认为您可以做出的最重要的性能改进是预先计算当前在 foreach 循环中重复计算的表达式。您可以使用单独的 foreach 循环来完成此操作。我还建议以不同的方式使用 solve 以避免第二次矩阵乘法:

X_gamma_list <- foreach(g=1:nrow(set)) %dopar% {
  X_gamma <- X[, which(set[g,] != 0)]
  I - (c/(1+c)) * (X_gamma %*% solve(crossprod(X_gamma), t(X_gamma)))
}

这些计算现在只执行一次，而不是对 Y 的每一列执行一次，这在您的情况下减少了 700 倍的工作量。

同样，按照 tim riffe 的建议，分解表达式 ((1+c)^(-sum(set[g,])/2)) 是有意义的，以及 -T/2 当我们这样做时:

a <- (1+c) ^ (-rowSums(set) / 2)
nT2 <- -T / 2

要迭代 zoo 对象 Y 的列，我建议使用 itertools 中的 isplitCols 函数> 包。确保在脚本顶部加载 itertools:

library(itertools)

isplitCols 让您仅发送每个任务所需的列，而不是将整个对象发送给所有工作人员。唯一的技巧是，您需要从生成的 zoo 对象中删除 dim 属性，代码才能正常工作，因为 isplitCols 使用 drop=TRUE。

最后，这是主 foreach 循环:

denom <- foreach(Yi=isplitCols(Y, chunkSize=1), .packages='zoo') %dopar% {
  dim(Yi) <- NULL  # isplitCols uses drop=FALSE
  result <- 0
  for(g in seq_along(X_gamma_list)) {
    S_g <- Yi %*% X_gamma_list[[g]] %*% Yi
    result <- result + a[g] * S_g ^ nT2
  }
  result
}

请注意，我不会并行执行内部循环。仅当 Y 中没有足够的列来保持所有处理器繁忙时，这才有意义。并行化内部循环可能会导致任务太短，从而有效地分块计算并使代码运行速度慢得多。由于 g 很大，因此有效地执行内部循环更为重要。