5

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

从这个公式中要注意的主要事情是S_g在矩阵乘法设置中使用 X_gamma 和 Y[,i]。X_gamma 取决于 value 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) 
 }

谢谢您的帮助！

4

2 回答 2

5

最明显的问题是您成为经典错误之一的受害者：没有预先分配输出向量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 倍的工作量。

类似地((1+c)^(-sum(set[g,])/2))，正如 tim riffe 所建议的那样，将表达式 分解出来是有意义的，以及-T / 2在我们处理它的时候：

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

要遍历zooobject的列Y，我建议使用包中的isplitCols函数itertools。itertools确保在脚本顶部加载：

library(itertools)

isplitCols让您只发送每个任务所需的列，而不是将整个对象发送给所有工作人员。唯一的技巧是您需要dim从结果对象中删除该属性zoo才能使您的代码正常工作，因为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。

于 2013-09-10T18:19:47.820 回答
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)为了速度而做。

于 2013-09-10T18:50:06.677 回答