R help - Aug 2013 - how to apply summation sign formula to R

Hi all

I have a short question relating to the usage of the summation sign in R.

Let's define A and B as two kxk matrice.
My goal is to calculate the matrix C for the periods from 1 to 200 (n=1-200).

C^(n) = ?_(j=1)^n [(?_(i=1)^(j-1) A^i ) B (?_(i=1)^(j-1) A^i)?  ]

How has that to be implemented in R (lets say for example for period = n = 150)?

Thank you for your help!!

Kind regards,
Sebastian

On Aug 24, 2013, at 6:37 AM, Sebastian Hersberger wrote:
> Hi all
> 
> I have a short question relating to the usage of the summation sign in R.
> 
> Let's define A and B as two kxk matrice.
> My goal is to calculate the matrix C for the periods from 1 to 200
(n=1-200).
> 
> C^(n) = ?_(j=1)^n [(?_(i=1)^(j-1) A^i ) B (?_(i=1)^(j-1) A^i)?  ]
> 
> How has that to be implemented in R (lets say for example for period = n =
150)?
I don't follow all this notation but you might want to look at the expm
package that has a powm function. It's possible that the expm function may
do what you want, but as I said I was unclear what the equation was attempting
to do. Perhaps if you explained the steps in mathematical language, then other
viewers who actually woke up to attend their linear algebra course might have a
better chance of offering assistance.

( When I saw the earlier copy yesterday, I had hopes that such wiser viewers
might chime in and still have such hopes.)

-- 

David Winsemius
Alameda, CA, USA

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