Folks:
1. Well, should have known better... I was wrong about what I said in
my earlier post below.
2. Here's yet another way to skin the cat using matrix multiplication
that is, maybe, a bit faster than using reduce.
First, the three methods (left as a homework exercise to see how/why they work):
## Originally proposed: d = dim of array test; xx is vector
e1 <- quote(Reduce('+', lapply(seq_len(d[3]), function(j) xx[j] *
test[, , j])))
## my earlier "better" (Ha!) way
e2 <- quote(apply(test*rep.int(xx, rep.int(d[1]*d[2],d[3])),1:2,sum))
## Using matrix multiplication
e3 <- quote(sapply(1:d[2],function(i)test[,i,]%*%xx))
Let's first try it with something simple to verify that the
expressions are "equivalent":
d <- 2:4
test <- array(1:24,dim=d)
xx <- 1:4
eval(e1)
eval(e2)
eval(e3)
IMPORTANT NOTE: This will **not** give identical() = TRUE results. e1
and e2 give integer results and e3 gives numeric. In general, floating
point arithmetic will defeat identical(). Instead, use all.equal()
Now here's a test bed.
d <- 10*c(50,50,4) ## dimensions of test array
test<- array(rnorm(prod(d)),dim=d)
xx <- runif(d[3])
results <- lapply(list(e1,e2,e3),function(e)system.time(eval(e)))
results
[[1]]
user system elapsed
0.34 0.01 0.36
[[2]]
user system elapsed
3.08 0.08 3.16
[[3]]
user system elapsed
0.27 0.00 0.26
##Try it with different dim
d <- 10*c(5,10,100) ## dimensions of test array
test<- array(rnorm(prod(d)),dim=d)
xx <- runif(d[3])
results <- lapply(list(e1,e2,e3),function(e)system.time(eval(e)))
results
[[1]]
user system elapsed
0.20 0.02 0.22
[[2]]
user system elapsed
0.25 0.01 0.27
[[3]]
user system elapsed
0.12 0.00 0.13
So I was dead wrong below; but matrix multiplication seems a tad
faster than Reduce.
Cheers to all,
-- Bert
On Mon, Mar 28, 2011 at 8:00 PM, Bert Gunter <bgun...@gene.com> wrote:
> Arrays are vectors stored in column major order. Hence
>
> H <- apply( H12*rep(hessian_lambda,e = n*k), 1:2, sum)
> ## but check that this gives what you want
>
> I'm pretty sure this is faster than the alternatives that were
> proposed, but whether it's much (or even noticeably) faster, I don't
> know.
>
>
>
> -- Bert
>
>
> On Mon, Mar 28, 2011 at 5:56 PM, Giuseppe <neo...@gmail.com> wrote:
>> Dennis,
>>
>> Yes, Reduce works, but unfortunately it does not speed up the
>> evaluation time (actually, it makes it a little bit slower).
>>
>> Thank you anyway.
>>
>>
>>
>> On Tue, Mar 29, 2011 at 2:21 AM, Dennis Murphy <djmu...@gmail.com> wrote:
>>> Hi:
>>>
>>> Is this what you're after?
>>>
>>> m <- array(1:27, c(3, 3, 3))
>>> xx <- 1:3
>>> Reduce('+', lapply(xx, function(k) xx[k] * m[, , k]))
>>>
>>>> m
>>> , , 1
>>>
>>> [,1] [,2] [,3]
>>> [1,] 1 4 7
>>> [2,] 2 5 8
>>> [3,] 3 6 9
>>>
>>> , , 2
>>>
>>> [,1] [,2] [,3]
>>> [1,] 10 13 16
>>> [2,] 11 14 17
>>> [3,] 12 15 18
>>>
>>> , , 3
>>>
>>> [,1] [,2] [,3]
>>> [1,] 19 22 25
>>> [2,] 20 23 26
>>> [3,] 21 24 27
>>>
>>>> xx <- 1:3
>>>> Reduce('+', lapply(xx, function(k) xx[k] * m[, , k]))
>>> [,1] [,2] [,3]
>>> [1,] 78 96 114
>>> [2,] 84 102 120
>>> [3,] 90 108 126
>>>
>>> HTH,
>>> Dennis
>>>
>>>
>>> On Mon, Mar 28, 2011 at 4:10 PM, Giuseppe <neo...@gmail.com> wrote:
>>>>
>>>> I have the following situation.
>>>>
>>>> H12 is an array of dimension (n,k,m) and hessian_lambda is a numeric
>>>> of length m.
>>>>
>>>> I need to multiply each matrix H12[,,1], H12[,,2], ..., H12[,m] by
>>>> hessian_lambda[1], hessian_lambda[2], ..., hessian_lambda[m],
>>>> respectively, and then add the resulting (n,k) matrices. What I am
>>>> using at the moment is the following code:
>>>>
>>>> H <- matrix(0, n, k)
>>>> for(j in 1:(m)) {
>>>> H <- H + H12[,,j]*hessian_lambda[j]
>>>> }
>>>>
>>>> Does anybody see a way to speed up this without using loop? (The
>>>> output needs to be a matrix).
>>>>
>>>> Thank you for your time.
>>>>
>>>> ______________________________________________
>>>> 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-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.
>>
>
>
>
> --
> Bert Gunter
> Genentech Nonclinical Biostatistics
> 467-7374
> http://devo.gene.com/groups/devo/depts/ncb/home.shtml
>
--
Bert Gunter
Genentech Nonclinical Biostatistics
467-7374
http://devo.gene.com/groups/devo/depts/ncb/home.shtml
______________________________________________
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.
