From: Prof Brian Ripley <ripley_at_stats.ox.ac.uk>

Date: Wed 05 Oct 2005 - 21:50:27 EST

Date: Wed 05 Oct 2005 - 21:50:27 EST

On Wed, 5 Oct 2005, Robin Hankin wrote:

> I have a square matrix Ainv of size N-by-N where N ~ 1000

*> I have a rectangular matrix H of size N by n where n ~ 4.
**> I have a vector d of length N.
**>
**> I need X = solve(t(H) %*% Ainv %*% H) %*% t(H) %*% Ainv %*% d
**>
**> and
**>
**> H %*% X.
**>
**> It is possible to rewrite X in the recommended crossprod way:
**>
**> X <- solve(quad.form(Ainv, H), crossprod(crossprod(Ainv, H), d))
**>
**> where quad.form() is a little function that evaluates a quadratic form:
**>
**> quad.form <- function (M, x){
**> jj <- crossprod(M, x)
**> return(drop(crossprod(jj, jj)))
**> }
*

That is not the same thing: it gives t(H) %*% Ainv %*% t(Ainv) %*% H .

**> QUESTION:
**

> 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.]

-- Brian D. Ripley, ripley@stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UK Fax: +44 1865 272595 ______________________________________________ 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 21:55:37 2005

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