Vectorizing many for loops

Question

I have a summation that I can calculate using four for loops but I wonder if this can be simplified, maybe using a vectorized function, to reduce the computation time. Something similar to the Kronecker product (in R: kronecker(x, x)), or maybe something using outer?

The summation is:

where E is the sample space of integers ranging from 1 - 9. The i and j indices are also integers ranging from 1 - 9. So, f, g, and h all matrices of dimension 9x9.

The h matrix is fixed and I have that but I am simulating g many times and then I choose the one that minimizes another function. The problem is, one thousand simulations, which is too few, takes about 1 second. I really want to try a million, but that many would take a long time.

I have the for loops in a function:

sim <- function(y, nreps, h) {

  G <- vector("list", nreps)     # list containing random values from Dirichlet distribution
  F <- vector("list", nreps)     # list containing the f matrices
  M <- vector("numeric", nreps)  # vector to store the results

  require(gtools)

  for(n in 1:nreps) {
    f <- matrix(0, nrow=9, ncol=9)        # initialize f
    g <- gtools::rdirichlet(9, rep(1,9))  # simulate g
    
    for(i in 1:9) {
      for(j in 1:9) {
        for(k in 1:9) {
          for(l in 1:9) {
            f[i,j] <- f[i,j] + h[i,k] * h[j,l] * g[k,l]  # summation (see above)
          }
        }
      }
    }
    F[[n]] <- f  # store f matrix
    G[[n]] <- g  # store g matrix
    M[n] <- sum((y - f)^2) # sum of squared differences between y and f
  }
  m <- which.min(M)  # which M is the minimum?
  return(list(g=G[[m]], m=M[m]))
}

And I call the function with

sim(y=f.y1, nreps=1000, h=x)

Here is the data:

> dput(f.y1)
structure(c(0.0182002022244692, 0.0121334681496461, 0.0101112234580384, 
0, 0, 0, 0, 0, 0, 0.0485338725985844, 0.0940343781597573, 0.112234580384226, 
0.0434782608695652, 0.00910010111223458, 0.00101112234580384, 
0, 0, 0, 0.0333670374115268, 0.110212335692619, 0.132457027300303, 
0.0808897876643074, 0.0222446916076845, 0.0070778564206269, 0.00101112234580384, 
0, 0, 0.0070778564206269, 0.0202224469160768, 0.0596562184024267, 
0.0616784630940344, 0.0262891809908999, 0.0070778564206269, 0, 
0, 0, 0.00202224469160768, 0.00505561172901921, 0.0151668351870576, 
0.0182002022244692, 0.0111223458038423, 0.00404448938321537, 
0, 0, 0, 0.00202224469160768, 0.00404448938321537, 0.00505561172901921, 
0.00505561172901921, 0.00202224469160768, 0.00202224469160768, 
0, 0, 0, 0, 0.00202224469160768, 0.00202224469160768, 0.00202224469160768, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0), class = "table", dim = c(9L, 9L), dimnames = structure(list(
    c("0", "1", "2", "3", "4", "5", "6", "7", "8"), c("0", "1", 
    "2", "3", "4", "5", "6", "7", "8")), names = c("", "")))

> dput(x)
structure(c(0.61, 0.16, 0.03, 0.005, 0, 0, 0, 0, 0, 0.32, 0.61, 
0.16, 0.03, 0.005, 0, 0, 0, 0, 0.06, 0.16, 0.61, 0.16, 0.03, 
0.005, 0, 0, 0, 0.01, 0.06, 0.16, 0.61, 0.16, 0.03, 0.01, 0, 
0, 0, 0.01, 0.03, 0.16, 0.61, 0.16, 0.03, 0.01, 0, 0, 0, 0.01, 
0.03, 0.16, 0.61, 0.16, 0.06, 0.01, 0, 0, 0, 0.005, 0.03, 0.16, 
0.61, 0.16, 0.06, 0, 0, 0, 0, 0.005, 0.03, 0.16, 0.61, 0.32, 
0, 0, 0, 0, 0, 0.005, 0.03, 0.16, 0.61), dim = c(9L, 9L))

And you'll need to load the gtools package for the rdirichlet function. Thanks heaps!

library(gtools)

Answer 1

Luckily this particular example is just "simple" matrix multiplication, so can easily be vectorised with:

sim1 <- function(y, nreps, h) {

  G <- vector("list", nreps)     # list containing random values from Dirichlet distribution
  F <- vector("list", nreps)     # list containing the f matrices
  M <- vector("numeric", nreps)  # vector to store the results

  require(gtools)

  for(n in 1:nreps) {
    g <- gtools::rdirichlet(9, rep(1,9))  # simulate g
    f <- h %*% g %*% t(h)
    
    F[[n]] <- f  # store f matrix
    G[[n]] <- g  # store g matrix
    M[n] <- sum((y - f)^2) # sum of squared differences between y and f
  }
  m <- which.min(M)  # which M is the minimum?
  return(list(g=G[[m]], m=M[m]))
}

Run function for comparison

#Original version
set.seed(0)
system.time(a <- sim(y=f.y1, nreps=1000, h=x))
#   user  system elapsed 
#   0.97    0.03    1.00 

#revised version
set.seed(0)
system.time(b <- sim1(y=f.y1, nreps=1000, h=x))
#   user  system elapsed 
#   0.01    0.00    0.02 

#Check they give the same answer
all.equal(a, b)  
#[1] TRUE