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