#setup options(width=60) olocale=Sys.setlocale(locale="C") library(Matrix) # sufficient statistics # y = Xb + e data(KNex) attach(KNex) ls(2) mmm <- as.matrix(mm) object.size(mmm)/1024 dim(mmm) length(y) # B(v) = (t(X) %*% X)^-1 t(X) %*% y(v) # # suff. statistics xx <- t(mmm) %*% mmm xy <- t(mmm) %*% y dim(xx) length(xy) b <- lm.fit(mmm,y)\$coef b.inv <- solve(t(mmm) %*% mmm) %*% t(mmm) %*% y all.equal(b,b.inv[,1],check.attr=FALSE,check.names=FALSE) # 30 permutations of y np <- 30 y.perm <- replicate(np,sample(y)) # use lm.fit system.time(for (i in 1:np) lm.fit(mmm,y.perm[,i])\$coef) # straight way system.time(for (i in 1:np) solve(t(mmm) %*% mmm) %*% t(mmm)%*%y.perm[,i]) # use solve, not inversion # solution above using inverse is one way of solving normal equations for b # (t(X) %*% X) %*% b = t(X) %*% y # there are other more efficient and stable ways, so let R solve it system.time(for (i in 1:np) solve(t(mmm) %*% mmm,t(mmm)%*%y.perm[,i])) # recycle suff. stat t(X) %*% X system.time({xx <- t(mmm) %*% mmm; for (i in 1:np) solve(xx,t(mmm) %*% y.perm[,i])}) # use crossprod functions instead system.time({xx <- crossprod(mmm); for (i in 1:np) solve(xx,crossprod(mmm,y.perm[,i]))}) # taking advantage of structure # xx = t(X) %*% X is positive definite (sometimes not, but leave that for another trickathon) # positive definite matrices have Cholesky factorizations # cholesky factor r satisfies # t(r) %*% r = t(mmm) %*% mmm # # but r is upper triangular # makes solving linear system much easier # to solve (t(r) %*% r) %*% b = z # solve t(r) %*% u = z (forward substitution) # solve r %*% b = u (backward substitution) # to solve using cholesky factorization r <- chol(crossprod(mmm)) u <- forwardsolve(r,crossprod(mmm,y),upper=TRUE,transpose=TRUE) b.chol <- backsolve(r,u) all.equal(b,b.chol[,1],check.attr=FALSE,check.names=FALSE) # what should really be recycled is the cholesky factor of sufficient statistic system.time({r <- chol(crossprod(mmm)); for (i in 1:np) backsolve(r, forwardsolve(r, crossprod(mmm,y.perm[,i]), upper=TRUE,trans=TRUE))}) # now using Matrix mm2 <- as(mm,"dgeMatrix") mm2[1:5,1:5] xx <- crossprod(mm2) class(xx) help("dpoMatrix-class") getClassDef("dpoMatrix") xx[1:10,1:10] mm[1:10,1:10] object.size(mm2)/1024 object.size(mm)/1024 b.Mat <- solve(crossprod(mm),crossprod(mm,y)) all.equal(b,b.Mat[,1],check.attr=FALSE,check.names=FALSE) # using solve system.time(for (i in 1:np) solve(crossprod(mm),crossprod(mm,y.perm[,i]))) # reusing suff. stat (reuses cholesky factor in the background) system.time({xx <- crossprod(mm); for (i in 1:np) solve(xx,crossprod(mm,y.perm[,i]))}) # the High Performance and Parallel Computing cran task view is here: # http://cran.r-project.org/web/views/HighPerformanceComputing.html