You may download the pdf Rmatlab
n=10000;
k=3;
e=rnorm(n);
x<-cbind(matrix(1,n,1),matrix(rnorm(n*k),n,k));
b=matrix(1,k+1,1);
y=x%*%b+e;
# using Choleskty solsution
xx<-crossprod(x,x)%*%solve(diag(k+1))
beta.chol<-crossprod(xx,crossprod(x,y))
system.time((beta.chol))
summary(lm(y~x))
