Dear All, I have a problem that I think could be solved much more efficiently, but I don't have a clue how to accomplish this. I have a matrix W with dimensions k*(p+1): Let's say W is 5*4 and looks like this:> W[,1] [,2] [,3] [,4] [1,] 1 -1 -1 1 [2,] 1 1 1 1 [3,] 1 2 -2 -1 [4,] 1 0 -1 -1 [5,] 1 -2 -1 0 I want to take each row (a p*1 vector), transpose it (1*p column vector), and then multiply that by the same row vector, which will give me a p*p matrix. So, for each row, I calculate: as.matrix(W[i,]) %*% W[i,] for i = 1, ..., k. Next, I want to find the sum of the k (p*p) matrices. The following code seems to work fine: k <- 5 p <- 3 W <- matrix(c( rep(1, k), round(rnorm(n=(p*k), mean=0, sd=1), 0)), ncol=(p+1), nrow=k) # create a random W matrix m1 <- matrix(0, nrow=p+1, ncol=p+1) for (i in 1:k) { m1 <- m1 + as.matrix(W[i,]) %*% W[i,] } but this just seems terribly inefficient. Any suggestions for improving this? Thanks in advance. Wolfgang -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
The "real" answer will tell us something about R syntax, arrays, vectorization, etc. I would also love to know it. But as it happens, the function you are describing is equal to t(W)%*%W. --- Wolfgang Viechtbauer <wviechtb at s.psych.uiuc.edu> wrote:> Dear All, > > I have a problem that I think could be solved much > more efficiently, but > I don't have a clue how to accomplish this. I have a > matrix W with > dimensions k*(p+1): > > Let's say W is 5*4 and looks like this: > > > W > [,1] [,2] [,3] [,4] > [1,] 1 -1 -1 1 > [2,] 1 1 1 1 > [3,] 1 2 -2 -1 > [4,] 1 0 -1 -1 > [5,] 1 -2 -1 0 > > I want to take each row (a p*1 vector), transpose it > (1*p column > vector), and then multiply that by the same row > vector, which will give > me a p*p matrix. So, for each row, I calculate: > > as.matrix(W[i,]) %*% W[i,] for i = 1, ..., k. > > Next, I want to find the sum of the k (p*p) > matrices. The following code > seems to work fine: > > > k <- 5 > p <- 3 > W <- matrix(c( rep(1, k), round(rnorm(n=(p*k), > mean=0, sd=1), 0)), > ncol=(p+1), nrow=k) # create a random W > matrix > m1 <- matrix(0, nrow=p+1, ncol=p+1) > for (i in 1:k) { > m1 <- m1 + as.matrix(W[i,]) %*% W[i,] > } > > > but this just seems terribly inefficient. Any > suggestions for improving > this? Thanks in advance. > > Wolfgang > > >-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-> r-help mailing list -- Read > http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html > Send "info", "help", or "[un]subscribe" > (in the "body", not the subject !) To: > r-help-request at stat.math.ethz.ch >_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._ __________________________________________________ -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Roger Peng suggests:> -----Original Message----- > From: Roger Peng [mailto:rpeng at stat.ucla.edu] > Sent: Saturday, August 10, 2002 2:26 PM > To: Wolfgang Viechtbauer > Cc: r-help at stat.math.ethz.ch > Subject: Re: [R] Help with improving efficiency > > In your example, you might consider > > m1 <- t(W) %*% diag(v) %*% W >[WNV] Mathematically correct but potentially memory demanding and slow. What's wrong with m1 <- crossprod(v*W, W) or, if you can assume that the elements of v are non-negative: m1 <- crossprod(sqrt(v) * W) ? I have learnt to be very wary about code that uses expliticly stored diagonal matrices. Bill Venables.> -roger > _______________________________ > UCLA Department of Statistics > rpeng at stat.ucla.edu > http://www.stat.ucla.edu/~rpeng > > On Fri, 9 Aug 2002, Wolfgang Viechtbauer wrote: > > > > The "real" answer will tell us something about R > > > syntax, arrays, vectorization, etc. I would also love > > > to know it. > > > > > > But as it happens, the function you are describing is > > > equal to t(W)%*%W. > > > > I must admit that I completely didn't notice that my code is just a > > cumbersome way of calculating the crossproduct of the W matrix (which > > could be accomplished even faster with the crossprod() function, as one > > person pointed out to me off-line). > > > > However, I guess I should have mentioned that the code within the loop > > is just a subset of the operations that are being performed on the row > > vectors. > > > > So, for example, let's say that each of the k (p*p) matrices is also > > being multiplied by a value v[i]: > > > > > > k <- 5 > > p <- 3 > > W <- matrix(c( rep(1, k), round(rnorm(n=(p*k), mean=0, sd=1), 0)), > > ncol=(p+1), nrow=k) # create a random W matrix > > v <- rnorm(n=k, mean=0, sd=1) # create random v values > > m1 <- matrix(0, nrow=p+1, ncol=p+1) > > for (i in 1:k) { > > m1 <- m1 + v[i] * as.matrix(W[i,]) %*% W[i,] > > } > > > > > > Obviously, this is different from simlpy taking the crossproduct of W. > > The best way to optimize the code probably depends on the types of > > operations being performed within the loop, but I am hoping to gain a > > general idea of how to optimize this particular code and then maybe I > > can apply that knowledge to the actual problem I am dealing with (and > > the "real" answer, as pointed out, might serve a nice illustration of > > how to write efficient code in general). > > > > Wolfgang > > > > > > > -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-. > -.-.- > > r-help mailing list -- Read > http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html > > Send "info", "help", or "[un]subscribe" > > (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch > > > _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._. > _._._ > > > > -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-. > -.-.- > r-help mailing list -- Read > http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html > Send "info", "help", or "[un]subscribe" > (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch > _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._. > _._._-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._