[R] how to apply summation sign formula to R

[Berend Hasselman]

On 25-08-2013, at 08:45, Prof Brian Ripley <ripley at stats.ox.ac.uk> wrote:

> On 25/08/2013 06:12, Berend Hasselman wrote:
>> 
>> On 24-08-2013, at 23:13, Sebastian Hersberger <sebastian.hersberger at unibas.ch> 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