Vectorized code for generating the Kac (Clement) matrix

On 26-04-2013, at 14:42, Enrico Schumann <es at enricoschumann.net> wrote:

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

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