R help - Oct 2005 - eliminate t() and %*% using crossprod() and solve(A,b)

Hi

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


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





--
Robin Hankin
Uncertainty Analyst
National Oceanography Centre, Southampton
European Way, Southampton SO14 3ZH, UK
  tel  023-8059-7743

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:
>
> 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)
> 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 at 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

you might want to see

Douglas Bates. Least squares calculations in R. R News,
4(1):17-20, June 2004.   http://CRAN.R-project.org/doc/Rnews/

he gives some rules of thumb, eg

use solve(A,b) not solve(A) %*% b
use crossprod(X) not t(X) %*% X
use crossprod(X,y) not t(X) y

----
Katharine Mullen
Department of Physics and Astronomy
Faculty of Sciences
Vrije Universiteit
de Boelelaan 1081
1081 HV Amsterdam
The Netherlands
room: T.1.06
tel: +31 205987870
fax: +31 205987992
e-mail: kate at nat.vu.nl
http://www.nat.vu.nl/~kate/