From: Robin Hankin <r.hankin_at_noc.soton.ac.uk>

Date: Wed 05 Oct 2005 - 23:08:05 EST

Date: Wed 05 Oct 2005 - 23:08:05 EST

On 5 Oct 2005, at 12:15, Dimitris Rizopoulos wrote:

*> Hi Robin,
**>
**> I've been playing with your question, but I think these two lines
**> are not equivalent:
**>
**> N <- 1000
**> n <- 4
**> Ainv <- matrix(rnorm(N * N), N, N)
**> H <- matrix(rnorm(N * n), N, n)
**> d <- rnorm(N)
*

> quad.form <- function (M, x){

*> jj <- crossprod(M, x)
**> return(drop(crossprod(jj, jj)))
**> }
**>
**>
**> X0 <- solve(t(H) %*% Ainv %*% H) %*% t(H) %*% Ainv %*% d
**> X1 <- solve(quad.form(Ainv, H), crossprod(crossprod(Ainv, H), d))
**> all.equal(X0, X1) # not TRUE
**>
**>
**> which is the one you want to compute?
**>
*

These are not exactly equivalent, but:

Ainv <- matrix(rnorm(1e6),1e3,1e3)

H <- matrix(rnorm(4000),ncol=4)

d <- rnorm(1000)

X0 <- solve(t(H) %*% Ainv %*% H) %*% t(H) %*% Ainv %*% d X1 <- solve(quad.form(Ainv, H), crossprod(crossprod(Ainv, H), d)) X0 - X1

[,1]

[1,] 4.884981e-15 [2,] 3.663736e-15 [3,] -5.107026e-15 [4,] 5.717649e-15

which is pretty close.

On 5 Oct 2005, at 12:50, Prof Brian Ripley wrote:

*>
***>> QUESTION:
**

I couldn't supply any performance data because I couldn't figure out the correct R commands to calculate H %*% X without using %*% or t()!

I was just wondering if there were a way to calculate

H %*% solve(t(H) %*% Ainv %*% H) %*% t(H) %*% Ainv %*% d

without using t() or %*%. And there doesn't seem to be (my original question didn't make it clear that I don't have X precalculated).

My take-home lesson from Brian Ripley is that H %*% X is fast --but this is only useful to me if one has X precalculated, and in general I don't. But this discussion suggests to me that it might be a good idea to change my routines and calculate X anyway.

thanks again Prof Ripley and Dimitris Rizopoulos

very best wishes

Robin

>> system.time(for(i in 1:100) t(H) %*% Ainv)

*>>
**> [1] 2.19 0.01 2.21 0.00 0.00
**>
**>> system.time(for(i in 1:100) crossprod(H, Ainv))
**>>
**> [1] 1.33 0.00 1.33 0.00 0.00
**>
**> so each is quite fast and the difference is not great. However,
**> that is
**> not comparing %*% with crossprod, but t & %*% with crossprod.
**>
**> I get
**>
**>
**>> system.time(for(i in 1:1000) H %*% X)
**>>
**> [1] 0.05 0.01 0.06 0.00 0.00
**>
**> which is hardly 'slow' (60 us for %*%), especially compared to
**> forming X
**> in
**>
**>
**>> system.time({X = solve(t(H) %*% Ainv %*% H) %*% t(H) %*% Ainv %*%
**>> d})
**>>
**> [1] 0.04 0.00 0.04 0.00 0.00
**>
**> I would probably have written
**>
**>
**>> system.time({X <- solve(crossprod(H, Ainv %*% H), crossprod
**>> (crossprod(Ainv, H), d))})
**>>
**> 1] 0.03 0.00 0.03 0.00 0.00
**>
**> which is faster and does give the same answer.
**>
**> [BTW, I used 2.2.0-beta which defaults to gcFirst=TRUE.]
**>
**> --
*

-- Robin Hankin Uncertainty Analyst National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743 ______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.htmlReceived on Wed Oct 05 23:12:30 2005

*
This archive was generated by hypermail 2.1.8
: Fri 03 Mar 2006 - 03:40:36 EST
*