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

2013-04-26 Thread Ravi Varadhan 
Thank you, Berend and Enrico, for looking into this.  I did not think of 
Enrico's clever use of cbind() to form the subsetting indices. 

Best,
Ravi

From: Berend Hasselman [b...@xs4all.nl]
Sent: Friday, April 26, 2013 10:08 AM
To: Enrico Schumann
Cc: Ravi Varadhan; 'r-help@r-project.org'
Subject: Re: [R] Vectorized code for generating the Kac (Clement) matrix

On 26-04-2013, at 14:42, Enrico Schumann  wrote:

> On Thu, 25 Apr 2013, Ravi Varadhan  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

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

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

2013-04-26 Thread Berend Hasselman 

On 26-04-2013, at 14:42, Enrico Schumann  wrote:

> On Thu, 25 Apr 2013, Ravi Varadhan  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

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

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

2013-04-26 Thread Enrico Schumann 
On Thu, 25 Apr 2013, Ravi Varadhan  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

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

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

2013-04-26 Thread Berend Hasselman 

On 25-04-2013, at 17:18, Ravi Varadhan  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

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

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

2013-04-26 Thread Ravi Varadhan 
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
rvarad...@jhmi.edu
410-502-2619


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.