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

Question

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

这个公式的主要注意事项是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。对于大型向量，一次附加一个值可能非常低效。

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

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 中没有足够的列来保持所有处理器忙碌时才有意义。并行化内部循环可能会导致任务太短，从而有效地unhunking计算并使代码运行得更慢。由于g 很大，因此有效地执行内部循环更为重要。

Answer 2

我第二个@eddi，您应该提供一些对象，以便我们可以运行代码。以下言论均基于凝视：

1) 你可以将S_g 保存在一个预先分配的向量中，并在循环之外执行最后一行（((1+c)^(-sum(set[g,])/2)) * ((S_g)^(-T/2))），因为rowSums(set) 会给你你需要的东西。这将删除一个使用 g

的索引实例

2) 索引会减慢您的速度。不要使用which()。逻辑向量工作得很好。

3) -T/2 很危险。它可以表示-0.5。如果这就是您想要的，那么只需 1/sqrt(S_g_vec) 以提高速度。