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

Berend Hasselman bhh at xs4all.nl
Fri Apr 26 13:12:22 CEST 2013

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

```