From: Dimitris Rizopoulos <dimitris.rizopoulos_at_med.kuleuven.be>

Date: Wed 05 Oct 2005 - 23:37:59 EST

Dimitris Rizopoulos

Ph.D. Student

Biostatistical Centre

School of Public Health

Catholic University of Leuven

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.html Received on Wed Oct 05 23:52:26 2005

Date: Wed 05 Oct 2005 - 23:37:59 EST

an alternative is to use X2 below, which seems to be a little bit faster:

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)))}

####################### #######################

invisible({gc(); gc(); gc()})

system.time(for(i in 1:200){

X0 <- solve(t(H) %*% Ainv %*% H) %*% t(H) %*% Ainv %*% d }, gcFirst = TRUE)

invisible({gc(); gc(); gc()})

system.time(for(i in 1:200){

X1 <- solve(quad.form(Ainv, H), crossprod(crossprod(Ainv, H), d)) }, gcFirst = TRUE)

invisible({gc(); gc(); gc()})

system.time(for(i in 1:200){

tH.Ainv <- crossprod(H, Ainv)

X2 <- solve(tH.Ainv %*% H) %*% colSums(t(tH.Ainv) * d)
}, gcFirst = TRUE)

all.equal(X0, X1)

all.equal(X0, X2)

I hope this helps.

Best,

Dimitris

Dimitris Rizopoulos

Ph.D. Student

Biostatistical Centre

School of Public Health

Catholic University of Leuven

Address: Kapucijnenvoer 35, Leuven, Belgium

Tel: +32/(0)16/336899 Fax: +32/(0)16/337015 Web: http://www.med.kuleuven.be/biostat/ http://www.student.kuleuven.be/~m0390867/dimitris.htm

- Original Message ----- From: "Robin Hankin" <r.hankin@noc.soton.ac.uk> To: "Prof Brian Ripley" <ripley@stats.ox.ac.uk> Cc: "RHelp" <r-help@stat.math.ethz.ch>; "Robin Hankin" <r.hankin@noc.soton.ac.uk> Sent: Wednesday, October 05, 2005 3:08 PM Subject: Re: [R] eliminate t() and %*% using crossprod() and solve(A,b)

*>
*

> 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:
**>>>
**>>> how to calculate
**>>>
**>>> H %*% X
**>>>
**>>> in the recommended crossprod way? (I don't want to take a
**>>> transpose
**>>> because t() is expensive, and I know that %*% is slow).
**>>>
**>>
**>> Have you some data to support your claims? Here I find (for random
**>> matrices of the dimensions given on a machine with a fast BLAS)
**>>
**>>
**>
**> 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.html
**>
*

Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm

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.html Received on Wed Oct 05 23:52:26 2005

*
This archive was generated by hypermail 2.1.8
: Sun 23 Oct 2005 - 18:20:41 EST
*