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

[Enrico Schumann]

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