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

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.

我的主要问题是IN REALITY 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倍.

在类似的情况下,((1+c)^(-sum(set[g,])/2))正如tim riffe所建议的,以及-T / 2在我们处理时,将表达式分解是有意义的:

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

为了迭代zoo对象的列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.