On 25-08-2013, at 08:45, Prof Brian Ripley <[email protected]> wrote:
> On 25/08/2013 06:12, Berend Hasselman wrote: >> >> On 24-08-2013, at 23:13, Sebastian Hersberger >> <[email protected]> wrote: >> >>> Thanks. I restate my problem/question and hope its better understandable >>> now. >>> >>> Let us define A and B as kxk matrices. C is the output (matrix), which I >>> try to calculate for differnt i values. >>> >>> So for example: I want to caluclate the matrix C for the value i=10: >>> >>> Therefore, I set: >>> >>> i <- c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10) >>> >>> Finally, I have to define the summation formula in R. My question is how >>> this following summation formula has to be applied to R. >>> >>> The arithmetic form of the formula equals: >>> >>> C = (Σ(from i=0 to i) A^i ) x B x (Σ(from i=0 to i) A^i )’ >>> >>> Which means: >>> matrix C equals the sum from i=0 to i times matrix A to the power of i >>> times matrix B >>> times the transposed/invers of the sum from i=0 to i times matrix A to the >>> power of i >> >> >> This is not the same (inner) summation as in the first post where i starts >> at 1 and goes to j-1. >> >> Original: (Σ_(i=1)^(j-1) A^i ) B (Σ_(i=1)^(j-1) A^i)’ >> That made me wonder what is supposed to happen when j=1? (Originally j >> started at 1 and stopped at n) >> >> David's solution can be wrapped in a function like this >> >> genAsum <- function(A,n,B) { >> tmp <- Reduce("+", lapply(0:n, FUN=function(j) A%^%j)) >> tmp %*% B %*% t(tmp) >> } > > It can, but as successive powers are best done by repeated multiplication (at > least for n as small as 10), a for() loop is preferable here. > > C <- diag(k); tmp <- matrix(0, k. k) > for(i in 0:n) { > tmp <- tmp + C > C <- C %*% A > } > > It is confused as to what i is here (it is used both as a dummy index and > apparently as what I have a 'n'). If you want this for n = 0, ..., 10 then > you should save intermediate results, and a for() loop is then even more > preferable. Your method is much faster even for large n. Example: set.seed(11) k <- 50 A <- matrix(runif(k*k),nrow=k) B <- matrix(runif(k*k),nrow=k) library(expm) g1 <- function(A,n,B) { tmp <- Reduce("+", lapply(0:n, FUN=function(j) A%^%j)) tmp %*% B %*% t(tmp) } n <- 100 D1 <- g1(A,n,B) mpow <- function(A,n) { k <- nrow(A) C <- diag(k); tmp <- matrix(0, k, k) for(i in 0:n) { tmp <- tmp + C C <- C %*% A } tmp } g2 <- function(A,n,B) { tmp <- mpow(A,n) tmp %*% B %*% t(tmp) } D2 <- g2(A,n,B) all.equal(D1,D2) # [1] TRUE library(compiler) g1.c <- cmpfun(g1) g2.c <- cmpfun(g2) library(rbenchmark) benchmark(g1(A,n,B),g2(A,n,B),g1.c(A,n,B),g2.c(A,n,B), replications=100, columns=c("test","replications","elapsed","relative")) # > benchmark(g1(A,n,B),g2(A,n,B),g1.c(A,n,B),g2.c(A,n,B), replications=100, # + columns=c("test","replications","elapsed","relative")) # test replications elapsed relative # 1 g1(A, n, B) 100 15.777 7.603 # 3 g1.c(A, n, B) 100 15.666 7.550 # 2 g2(A, n, B) 100 2.184 1.053 # 4 g2.c(A, n, B) 100 2.075 1.000 Berend ______________________________________________ [email protected] 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.
