We seek to compare different methods of fitting a linear regression.

First we create a matrix and vector to serve as our design matrix and response variables.

x <- matrix(0,256000,512)
      x[,1] <- 1
      for (l in 2:512) {
        x[,l] <- rnorm(256000)
        gc()
      }
# vector of parameters
b <- seq(1,512)

y <- x %*% b+rnorm(256000)

lm() function

The default method for fitting a linear model to a set of data. Uses “qr” method for fitting. This method may refer to the QR-decomposition method of solving for least squares, but needs a refernce to confirm. Calls on the function lm.fit() for numerical regression fitting.

system.time(lm(y~x))
##    user  system elapsed 
##   75.61    1.63   77.24

naive matrix operations

Using matrix multiplication for solving for least squares. \[B=(X'X)^{-1}X'Y\] Computes using elementary matrix functions, t() transpose, solve() inverse calculation, and %*% for matrix multiplication.

system.time(solve(t(x) %*% x) %*% t(x) %*% y)
##    user  system elapsed 
##   75.60    0.63   76.36

crossprod()

Uses matrix operations to solve least squares but calls upon the crossprod() function. crossprod(x,y) is the equivalent of t(x)%*%y. This function is supposedly quicker than the explicit transpose and multiplication operations.

system.time(solve(crossprod(x)) %*% crossprod(x,y))
##    user  system elapsed 
##   35.59    0.00   35.61

ginv()

Penrose-Moore Generalized Inverse. A method of finding a pseudo-inverse of a matrix whose properties make it suitable for solving least squares. Let \(X\) be a real matrix with linearly independent columns, then the Penrose-Moore pseudo-inverse \(X^+\) has the property, \[X^+=(X'X)^{-1}X'\] Thus \[B=X^+Y\] ginv() calls upon the singular value decomposition function, svd(), to find the pseudo-inverse. Let \(X=U\Sigma V^*\) be the svd factorization of \(X\), where \(^*\) is the complex conjugate operator. Then the pseudo-inverse is \[X^+=V\Sigma^+U^*\] Where \(\Sigma^+\) is the pseudo-inverse of \(\Sigma\) formed by replacing every non-zero diagonal entry by its reciprocal and transposing the resulting matrix. We can possibly find a more efficient means by using svd().

# display source code
library(MASS)
system.time(ginv(x) %*% y)
##    user  system elapsed 
##  219.98    1.36  221.77

qr decompostion

Orthogonal matrices preserve Euclidean distances, thus they are a good candidate for least squares solutions. Let \(X\) be a real matrix, then \(X=QR\) is the QR factorization of \(X\). From this form, the least squares solution is \[B=R^{-1}Q'Y\]

system.time(qr.solve(x,y))
##    user  system elapsed 
##   64.94    0.67   65.68

svd decomposition

In the same way that the ginv() function operates, we use the svd() function to find the least squares solution. \[B=V\Sigma^+U^*Y\]

Xsvd <- svd(x)
system.time(Xsvd$v %*% (1/Xsvd$d * t(Xsvd$u)))
##    user  system elapsed 
##   32.36    0.48   32.87
rm(Xsvd)

lm.fit()

lm.fit() is the lower level function that performs the numerical calculations in fitting regression. This function is called on by the lm() function. The default method for this function is set to “qr”.

system.time(lm.fit(x,y))
##    user  system elapsed 
##   65.63    0.19   65.86

lsfit()

Finds the least squares estimate \(B\) in the following model \[Y=XB+e\] Note there is not much information in the help page on this function. Method is most likely the QR-decomposition method.

system.time(lsfit(x,y,intercept = FALSE))
##    user  system elapsed 
##   66.37    0.50   67.00

.lm.fit()

“.lm.fit() is bare bone wrapper to the innermost QR-based C code”

system.time(.lm.fit(x,y))
##    user  system elapsed 
##   65.19    0.14   65.42