Ordinary Least Squares

• we want to solve something of the form \[ \boldsymbol{Y} = \boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{\epsilon} \]

Solutions

• OLS is typically learned as

\[ (\boldsymbol{X}'\boldsymbol{X})^{-1}\boldsymbol{X}'\boldsymbol{Y} \] - but in practice, this is computationally inefficient - floating point algebra, matrix inversion, etc etc - instead, we use QR decomposition

QR Decomposition

\[ \boldsymbol{X} = \boldsymbol{QR} \] - since R is an upper triangular matrix, we can set

\[ \boldsymbol{Y} = \boldsymbol{QR}\boldsymbol{\beta} \] - meaning that:

\[ \boldsymbol{Q}^{-1}\boldsymbol{Y} = \boldsymbol{R}\boldsymbol{\beta} \] - by definition of the QR decomposition: \[ \boldsymbol{Q}^{-1}\boldsymbol{Y} = \boldsymbol{Q}'\boldsymbol{Y}\] - we can therefore exploit the triangular structur of \(\boldsymbol{R}\) to backsolve the equation: \[ \boldsymbol{Q}'\boldsymbol{Y} = \boldsymbol{R}\boldsymbol{\beta} \]

LU Decomposition

\[ \left(\boldsymbol{X}'\boldsymbol{X}\right)\boldsymbol{b} = \boldsymbol{X}'\boldsymbol{y} \]

Implementation

ols <- function(data, formula, method="QR") {
  
  mf <- model.frame(formula, data, na.action=na.exclude)
  Y <- model.response(mf, type="double")
  X <- model.matrix(formula, mf)
  
  if (method == "matrix") b <- solve(crossprod(X)) %*% crossprod(X, Y)
  if (method == "QR") {
    QR <- qr(X)
    Q <- qr.Q(QR)
    R <- qr.R(QR)
    b <- backsolve(R, crossprod(Q,y))
    rownames(b) <- colnames(X)
  }
  
  return(b)
}

Testing the function:

set.seed(42)

N <- 1000
X1 <- rnorm(N)
X2 <- rnorm(N)
y = 1.5 + 2*X1 + 3*X2 + rnorm(N)

df <- data.frame(
  y, X1, X2
)

f <- y ~ X1 + X2

coef(lm(f, df))

## (Intercept)          X1          X2 
##    1.496133    1.978201    2.981712

ols(df, f)

##                 [,1]
## (Intercept) 1.496133
## X1          1.978201
## X2          2.981712

ols(df, f, "matrix")

##                 [,1]
## (Intercept) 1.496133
## X1          1.978201
## X2          2.981712