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