PROBLEM SET 1

Consider the unsolvable \(Ax = b\) as given below:

\[ \left[ {\begin{array}{cc} 1 & 0\\ 1 & 1\\ 1 & 3\\ 1 & 4\\ \end{array} } \right] \left[ {\begin{array}{c} x_1\\ x_2\\ \end{array} } \right] = \left[ {\begin{array}{c} 0\\ 8\\ 8\\ 19\\ \end{array} } \right]\]

Write code in R to compute \(A^T A\) and \(A^T b\).

A<-matrix(c(1,1,1,1,0,1,3,4),nrow=4)
b<- matrix(c(0,8,8,20), nrow=4)
A_T_A <- t(A) %*% A
A_T_b <- t(A) %*% b
A_T_A 
##      [,1] [,2]
## [1,]    4    8
## [2,]    8   26
A_T_b
##      [,1]
## [1,]   36
## [2,]  112

Solve for \(\hat{x}\) in R using the above two computed matrices.

x_hat <- solve(A_T_A) %*% A_T_b
x_hat
##      [,1]
## [1,]    1
## [2,]    4
#sanity test, removing 0
lsfit(A, b,intercept=FALSE)$coef
## X1 X2 
##  1  4

What is the squared error of this solution?

p <- A %*% x_hat
#error vector
e <- b-p
#sanity test
e
##      [,1]
## [1,]   -1
## [2,]    3
## [3,]   -5
## [4,]    3
lsfit(A, b,intercept=FALSE)$residuals
## [1] -1  3 -5  3
# Find the squared error
error_sq <- sum(e^2)
error_sq
## [1] 44

Instead of b = [0;8;8;20], start with p = [1;5;13;17] and find the exact solution.

p <- matrix(c(1,5,13,17), nrow=4) 
A_T_p <- t(A) %*% p
x_hat2 <- solve(A_T_A) %*% A_T_p
p2 <- A %*% x_hat2
p2
##      [,1]
## [1,]    1
## [2,]    5
## [3,]   13
## [4,]   17
e2 <- p2-p
e2
##              [,1]
## [1,] 0.000000e+00
## [2,] 8.881784e-16
## [3,] 3.552714e-15
## [4,] 3.552714e-15
#We see that e2 is essentially 0
error_sq2 <- round(sum(e2^2))
error_sq2 
## [1] 0

Show that the error \(e = b - p = [-1; 3;-5; 3]\)

#residual error
b-p
##      [,1]
## [1,]   -1
## [2,]    3
## [3,]   -5
## [4,]    3
#sanity test
lsfit(A, b,intercept=FALSE)$residuals
## [1] -1  3 -5  3

Show that the error \(e\) is orthogonal to \(p\) and to each of the columns of A.

# e*p == 0
sum(e*p)
## [1] -1.030287e-13
#e and the columns of A are orthogonal
#column 1
sum(e*A[,1])
## [1] -7.993606e-15
#column 2
sum(e*A[,2])
## [1] -2.575717e-14

PROBLEM SET 2

Consider the modified auto-mpg data (obtained from the UC Irvine Machine Learning dataset). This dataset contains 5 columns: displacement, horsepower, weight, acceleration, mpg. We are going to model mpg as a function of the other four variables.

#set working directory to current
setwd_thisdir <- function () {
  this.dir <- dirname(parent.frame(3)$ofile)
  setwd(this.dir)
}
#read data into matrix format
autoData<-as.matrix(read.table("auto-mpg.data", header=FALSE,as.is=TRUE))
colnames(autoData) <- c("displacement", "horsepower","weight","acceleration","mpg")

Extract an A matrix from the first four columns and vector b from the fifth.

A <- data.matrix(autoData[,1:4])
b <- data.matrix(autoData[,5])
# Adding a column of ones for the Intercepts
intercept <- rep(1, nrow(A))
A <- cbind(intercept, A)

Using the least squares approach, find the best fitting solution:

#using approach from Problem Set 1
A_T_A <- t(A) %*% A
A_T_b <- t(A) %*% b
#solve for x - hat
x_hat <- solve(A_T_A) %*% A_T_b
#or
#lsfit(A, b,intercept=FALSE)$coef
#squared error
p <- A %*% x_hat
#error vector
e <- b-p
#or
#lsfit(A, b,intercept=FALSE)$residuals
# Find the squared error
error_sq <- sum(e^2)
error_sq
## [1] 6979.413
#check results against xhat
lsfit(A, b, intercept=FALSE)$coefficients
##    intercept displacement   horsepower       weight acceleration 
## 45.251139699 -0.006000871 -0.043607731 -0.005280508 -0.023147999
x_hat
##                      [,1]
## intercept    45.251139699
## displacement -0.006000871
## horsepower   -0.043607731
## weight       -0.005280508
## acceleration -0.023147999