[R] Matrix oriented computing

[Gabor Grothendieck]

On 8/26/05, Marc Schwartz (via MN) <mschwartz at mn.rr.com> wrote:
> On Fri, 2005-08-26 at 15:25 +0200, Peter Dalgaard wrote:
> > Marc Schwartz <MSchwartz at mn.rr.com> writes:
> >
> > > x <- c(0.005, 0.010, 0.025, 0.05, 0.1, 0.5, 0.9,
> > >        0.95, 0.975, 0.99, 0.995)
> > >
> > > df <- c(1:100)
> > >
> > > mat <- sapply(x, qchisq, df)
> > >
> > > > dim(mat)
> > > [1] 100  11
> > >
> > > > str(mat)
> > >  num [1:100, 1:11] 3.93e-05 1.00e-02 7.17e-02 2.07e-01 4.12e-01 ...
> >
> > outer() is perhaps a more natural first try... It does give the
> > transpose of the sapply approach, though.
> >
> > round(t(outer(x,df,qchisq)),2)
> >
> > should be close. You should likely add dimnames.
> 
> 
> 
> What I find interesting, is that I would have intuitively expected
> outer() to be faster than sapply().  However:
> 
> 
> >  system.time(mat <- sapply(x, qchisq, df), gcFirst = TRUE)
> [1] 0.01 0.00 0.01 0.00 0.00
> 
> >  system.time(mat1 <- round(t(outer(x, df, qchisq)), 2),
>               gcFirst = TRUE)
> [1] 0.01 0.00 0.01 0.00 0.00
> 
> # No round() or t() to test for overhead
> >  system.time(mat2 <- outer(x, df, qchisq), gcFirst = TRUE)
> [1] 0.01 0.00 0.02 0.00 0.00
> 
> 
> # Bear in mind the round() on mat1 above
> > all.equal(mat, mat1)
> [1] "Mean relative  difference: 4.905485e-05"
> 
> > all.equal(mat, t(mat2))
> [1] TRUE
> 
> 
> Even when increasing the size of 'df' to 1:1000:
> 
> 
> >  system.time(mat <- sapply(x, qchisq, df), gcFirst = TRUE)
> [1] 0.16 0.01 0.16 0.00 0.00
> 
> >  system.time(mat1 <- round(t(outer(x, df, qchisq)), 2), gcFirst =
> TRUE)
> [1] 0.16 0.00 0.18 0.00 0.00
> 
> >  # No round() or t() to test for overhead
> >  system.time(mat2 <- outer(x, df, qchisq), gcFirst = TRUE)
> [1] 0.16 0.01 0.17 0.00 0.00
> 
> 
> 
> It also seems that, at least in this case, t() and round() do not add
> much overhead.
> 

You might need to do it repeatedly to get a more reliable reading.
When I do it 1000 times outer is faster than sapply though not by much:

> n <- 1000
> system.time(for (i in 1:n) mat <- sapply(x, qchisq, df), gcFirst = TRUE)
[1] 14.05  0.00 14.43    NA    NA
> 
> system.time(for(i in 1:n) mat2 <- outer(x, df, qchisq), gcFirst = TRUE)
[1] 13.42  0.00 13.85    NA    NA