[R] eliminate t() and %*% using crossprod() and solve(A,b)

[Robin Hankin]

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