Thank you, Berend and Enrico, for looking into this. I did not think of Enrico's clever use of cbind() to form the subsetting indices.
Best, Ravi ________________________________________ From: Berend Hasselman [b...@xs4all.nl] Sent: Friday, April 26, 2013 10:08 AM To: Enrico Schumann Cc: Ravi Varadhan; 'r-help@r-project.org' Subject: Re: [R] Vectorized code for generating the Kac (Clement) matrix On 26-04-2013, at 14:42, Enrico Schumann <e...@enricoschumann.net> wrote: > On Thu, 25 Apr 2013, Ravi Varadhan <ravi.varad...@jhu.edu> writes: > >> Hi, I am generating large Kac matrices (also known as Clement matrix). >> This a tridiagonal matrix. I was wondering whether there is a >> vectorized solution that avoids the `for' loops to the following code: >> >> n <- 1000 >> >> Kacmat <- matrix(0, n+1, n+1) >> >> for (i in 1:n) Kacmat[i, i+1] <- n - i + 1 >> >> for (i in 2:(n+1)) Kacmat[i, i-1] <- i-1 >> >> The above code is fast, but I am curious about vectorized ways to do this. >> >> Thanks in advance. >> Best, >> Ravi >> > > > This may be a bit faster; but as Berend and you said, the original > function seems already fast. > > n <- 5000 > > f1 <- function(n) { > Kacmat <- matrix(0, n+1, n+1) > for (i in 1:n) Kacmat[i, i+1] <- n - i + 1 > for (i in 2:(n+1)) Kacmat[i, i-1] <- i-1 > Kacmat > } > f3 <- function(n) { > n1 <- n + 1L > res <- numeric(n1 * n1) > dim(res) <- c(n1, n1) > bw <- n:1L ## bw = backward, fw = forward > fw <- seq_len(n) > res[cbind(fw, fw + 1L)] <- bw > res[cbind(fw + 1L, fw)] <- fw > res > } > > system.time(K1 <- f1(n)) > system.time(K3 <- f3(n)) > identical(K3, K1) > > ## user system elapsed > ## 0.132 0.028 0.161 > ## > ## user system elapsed > ## 0.024 0.048 0.071 > ## Using some of your code in my function I was able to speed up my function f2. Complete code: f1 <- function(n) { #Ravi Kacmat <- matrix(0, n+1, n+1) for (i in 1:n) Kacmat[i, i+1] <- n - i + 1 for (i in 1:n) Kacmat[i+1, i] <- i Kacmat } f2 <- function(n) { # Berend 1 modified to use 1L Kacmat <- matrix(0, n+1, n+1) Kacmat[row(Kacmat)==col(Kacmat)-1L] <- n:1L Kacmat[row(Kacmat)==col(Kacmat)+1L] <- 1L:n Kacmat } f3 <- function(n) { # Enrico n1 <- n + 1L res <- numeric(n1 * n1) dim(res) <- c(n1, n1) bw <- n:1L ## bw = backward, fw = forward fw <- seq_len(n) res[cbind(fw, fw + 1L)] <- bw res[cbind(fw + 1L, fw)] <- fw res } f4 <- function(n) {# Berend 2 using which with arr.ind=TRUE Kacmat <- matrix(0, n+1, n+1) k1 <- which(row(Kacmat)==col(Kacmat)-1L, arr.ind=TRUE) k2 <- which(row(Kacmat)==col(Kacmat)+1L, arr.ind=TRUE) Kacmat[k1] <- n:1L Kacmat[k2] <- 1L:n Kacmat } library(compiler) f1.c <- cmpfun(f1) f2.c <- cmpfun(f2) f3.c <- cmpfun(f3) f4.c <- cmpfun(f4) f1(n) f2(n) n <- 5000 system.time(K1 <- f1(n)) system.time(K2 <- f2(n)) system.time(K3 <- f3(n)) system.time(K4 <- f4(n)) system.time(K1c <- f1.c(n)) system.time(K2c <- f2.c(n)) system.time(K3c <- f3.c(n)) system.time(K4c <- f4.c(n)) identical(K2,K1) identical(K3,K1) identical(K4,K1) identical(K1c,K1) identical(K2c,K2) identical(K3c,K3) identical(K4c,K4) Result: # > system.time(K1 <- f1(n)) # user system elapsed # 0.387 0.120 0.511 # > system.time(K2 <- f2(n)) # user system elapsed # 3.541 0.702 4.250 # > system.time(K3 <- f3(n)) # user system elapsed # 0.108 0.089 0.199 # > system.time(K4 <- f4(n)) # user system elapsed # 1.975 0.355 2.336 # > # > system.time(K1c <- f1.c(n)) # user system elapsed # 0.323 0.120 0.445 # > system.time(K2c <- f2.c(n)) # user system elapsed # 3.374 0.422 3.807 # > system.time(K3c <- f3.c(n)) # user system elapsed # 0.107 0.098 0.205 # > system.time(K4c <- f4.c(n)) # user system elapsed # 1.816 0.384 2.203 # > identical(K2,K1) # [1] TRUE # > identical(K3,K1) # [1] TRUE # > identical(K4,K1) # [1] TRUE # > identical(K1c,K1) # [1] TRUE # > identical(K2c,K2) # [1] TRUE # > identical(K3c,K3) # [1] TRUE # > identical(K4c,K4) # [1] TRUE So Ravi's original and Enrico's versions are the quickest. Using which with arr.ind made my version run a lot quicker. All in all an interesting exercise. Berend ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.