R help - Apr 2013 - Vectorized code for generating the Kac (Clement) matrix

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

Ravi Varadhan, Ph.D.
Assistant Professor
The Center on Aging and Health
Division of Geriatric Medicine & Gerontology
Johns Hopkins University
rvaradhan@jhmi.edu<mailto:rvaradhan@jhmi.edu>
410-502-2619


	[[alternative HTML version deleted]]

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

On Thu, 25 Apr 2013, Ravi Varadhan <ravi.varadhan at 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 
## 


-- 
Enrico Schumann
Lucerne, Switzerland
http://enricoschumann.net