Effective Multiple Matrix Multiplication Method

Question

I am trying to do matrix multiplication S_g for every i, and every g with i. This is what I have tried so far, but it takes a huge amount of time to complete. Is there a more efficient computational method to do the same?

The main thing that follows from this formula: S_g uses X_gamma and Y [, i] in setting up matrix multiplication. X_gamma depends on the value of g . Therefore, for each I, I must perform g matrix multiplications.

Here is the logic:

• For each I, the calculation should be performed for each g. Then, for each g, X_gamma is selected as a subset of X. This is how X_gamma is determined. Take g = 3. When we look at "set [3,]", we get that column B is the only one with value! = 0. Therefore, I select column B in X, and that will be X_gamma.

My main problem: IN REALITY, g = 13,000 and 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) } 

Thank you for your help!

Answer 1

The most obvious problem is that you are the victim of one of the classic mistakes: do not prevail the output vector result . Adding one value at a time can be very inefficient for large vectors.

In your case, the result does not have to be a vector: you can accumulate the results in one value:

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

But I think that the most important performance improvement you could make is to precompile the expressions that are currently being re-evaluated in the foreach . You can do this with a separate foreach . I also suggest using solve in different ways to avoid a second matrix multiplication:

 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))) } 

These calculations are now performed only once, and not once for each Y column, which is 700 times less in your case.

In the same spirit, it makes sense to expand the expression ((1+c)^(-sum(set[g,])/2)) , as suggested by tim riffe, as well as -T / 2 , while we are in it:

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

To isplitCols over the columns of the zoo Y object, I suggest using the isplitCols function from the itertools package. Make sure you download itertools at the top of the script:

 library(itertools) 

isplitCols allow you to send only those columns that are necessary for each task, and not send the entire object to all employees. The only trick is that you need to remove the dim attribute from the resulting zoo objects for your code to work, since isplitCols uses drop=TRUE .

Finally, here is the main 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 } 

Please note that I will not execute the inner loop in parallel. This would only make sense if there weren’t enough columns in Y to keep all your processors busy. Parallelizing the inner loop can lead to tasks that are too short, effectively turn off computation, and make code much slower. It is much more important to efficiently execute the inner loop, since g large.

Answer 2

I am the second @eddi that you have to give some objects so that we can run the code. The following comments are based on the look:

1) you can save S_g in a pre-distributed vector and make this last line (( ((1+c)^(-sum(set[g,])/2)) * ((S_g)^(-T/2)) ) outside the loop since rowSums(set) will provide you with what you need. This removes one indexing instance using g

2) indexing slows you down. Do not use which() . Logical vectors work fine.

3) -T/2 is dangerous. This could mean -0.5 . If this is what you want, then just 1/sqrt(S_g_vec) for speed.