[R] Vectorized code for generating the Kac (Clement) matrix

[Berend Hasselman]

On 25-04-2013, at 17:18, Ravi Varadhan <ravi.varadhan at jhu.edu> wrote:

> 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.

You could vectorize like this

Kacmat <- matrix(0, n+1, n+1)  
Kacmat[row(Kacmat)==col(Kacmat)-1] <- n -(1:n) + 1
Kacmat[row(Kacmat)==col(Kacmat)+1] <- 1:n

But this show that your version is pretty quick

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
}

f2 <- function(n) {
    Kacmat <- matrix(0, n+1, n+1)  
    Kacmat[row(Kacmat)==col(Kacmat)-1] <- n -(1:n) + 1
    Kacmat[row(Kacmat)==col(Kacmat)+1] <-1:n
    Kacmat
}

library(compiler)

f1.c <- cmpfun(f1)
f2.c <- cmpfun(f2)

n <- 5000

system.time(K1 <- f1(n))        
system.time(K2 <- f2(n))        
system.time(K3 <- f1.c(n))        
system.time(K4 <- f2.c(n))        
identical(K2,K1)
identical(K3,K1)
identical(K4,K1)  

# > system.time(K1 <- f1(n))        
#    user  system elapsed 
#   0.386   0.120   0.512 
# > system.time(K2 <- f2(n))        
#    user  system elapsed 
#   3.779   1.141   4.940 
# > system.time(K3 <- f1.c(n))        
#    user  system elapsed 
#   0.323   0.119   0.444 
# > system.time(K4 <- f2.c(n))        
#    user  system elapsed 
#   3.607   0.852   4.472 
# > identical(K2,K1)
# [1] TRUE
# > identical(K3,K1)
# [1] TRUE
# > identical(K4,K1)
# [1] TRUE


Berend