Problem Set 1

In this problem set we’ll work out some properties of the least squares solution that we reviewed in the weekly readings. Consider the unsolvable system Ax = b as given below:

\[ \left[ \begin{array}{cc} 1 & 0 \\ 1 & 1 \\ 1 & 3 \\ 1 & 4 \end{array} \right] % \left[ \begin{array}{cc} x1 \\ x2 \end{array} \right] = \left[ \begin{array}{cc} 0 \\ 8 \\ 8 \\ 20 \end{array} \right] \]

Write R markdown script to compute \(A^TA\) and \(A^Tb\).

A <- matrix(c(1, 1, 1, 1, 0, 1, 3, 4), nrow = 4)
b <- c(0, 8, 8, 20)

# compute A^TA and A^Tb
(ata <- t(A) %*% A)
##      [,1] [,2]
## [1,]    4    8
## [2,]    8   26
(atb <- t(A) %*% b)
##      [,1]
## [1,]   36
## [2,]  112

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

As \(A^TA\) \(\hat{x}\) = \(A^Tb\), so \(\hat{x}\) = \(A^Tb\) \(A^TA^{-1}\)

# calculate xhat called leastSqr
(leastSqr <- solve(ata) %*% atb)
##      [,1]
## [1,]    1
## [2,]    4

What is the squared error of this solution?

\[b = \left[ \begin{array}{cc} 0 \\ 8 \\ 8 \\ 20 \end{array} \right] \]

and

\[ A \hat{x} = \left[ \begin{array}{cc} 1 & 0 \\ 1 & 1 \\ 1 & 3 \\ 1 & 4 \end{array} \right] % \left[ \begin{array}{cc} 1 \\ 4 \end{array} \right] = \left[ \begin{array}{cc} 1 \\ 5 \\ 13 \\ 7 \end{array} \right] \]

Hence errSqr = sum((b - A \(\hat{x})^2\))

# squred error
(errSqr <- sum((b - A %*% leastSqr)^2))
## [1] 44

Instead of b = [0; 8; 8; 20], start with p = [1; 5; 13; 17] and find the exact solution (i.e. show that this system is solvable as all equations are consistent with each other. This should result in an error vector e = 0).

p <- c(1, 5, 13, 17)
atp <- t(A) %*% p
(leastSqr <- solve(ata) %*% atp)
##      [,1]
## [1,]    1
## [2,]    4
# squred error
(errSqr <- sum((p - A %*% leastSqr)^2))
## [1] 2.603241e-29

Show that the error e = b - p = [-1; 3; -5; 3].

e <- c(-1, 3, -5, 3)
isTRUE(all.equal(e, (b - p)))
## [1] TRUE

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

# the error e is orthogonal to p and to each of the columns of A when the transpose of e 
# 0 means the angle between them is 90 degrees
t(e) %*% p
##      [,1]
## [1,]    0
t(e) %*% A[, 1]
##      [,1]
## [1,]    0
t(e) %*% A[, 2]
##      [,1]
## [1,]    0

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.

Write an R markdown script that takes in the auto-mpg data, extracts an A matrix from the first 4 columns and b vector from the fifth (mpg) column. Using the least squares approach, your code should compute the best fitting solution. That is, find the best fitting equation that expresses mpg in terms of the other 4 variables. Finally, calculate the fitting error between the predicted mpg of your model and the actual mpg. Your script should be able to load in the 5 column data set, extract A and b, and perform the rest of the calculations. Please have adequate comments in your code to make it easy to follow your work.

# read in data table, and rename
auto.mpg <- read.table('auto-mpg.data', quote='\'', comment.char='')
names(auto.mpg) <- c('displacement', 'horsepower', 'weight', 'acceleration', 'mpg')
head(auto.mpg)
##   displacement horsepower weight acceleration mpg
## 1          307        130   3504         12.0  18
## 2          350        165   3693         11.5  15
## 3          318        150   3436         11.0  18
## 4          304        150   3433         12.0  16
## 5          302        140   3449         10.5  17
## 6          429        198   4341         10.0  15
# create A and b
A <- data.matrix(auto.mpg[1:4])
b <- data.matrix(auto.mpg[5])
A <- cbind(1, A)

# compute A^TA and A^Tb
(ata <- t(A) %*% A)
##                        displacement  horsepower     weight acceleration
##                  392.0      76209.5     40952.0    1167213       6092.2
## displacement   76209.5   19097634.2   9374647.0  259345480    1123011.9
## horsepower     40952.0    9374647.0   4857524.0  132989885     607832.3
## weight       1167213.0  259345480.0 132989885.0 3757575489   17758103.6
## acceleration    6092.2    1123011.9    607832.3   17758104      97656.9
(atb <- t(A) %*% b)
##                     mpg
##                  9190.8
## displacement  1529685.9
## horsepower     868718.8
## weight       25209061.4
## acceleration   146401.4
# calculate xhat
(leastSqr <- solve(ata) %*% atb)
##                       mpg
##              45.251139699
## displacement -0.006000871
## horsepower   -0.043607731
## weight       -0.005280508
## acceleration -0.023147999
# calculate the squared error
(errSqr <- sum((b - A %*% leastSqr)^2))
## [1] 6979.413
# the least square equation for 'displacement', 'horsepower', 'weight', 'acceleration' expressing 'mpg'
paste0('mpg = ', round(leastSqr[2, 1], digits=2), '*displacement + (', 
                     round(leastSqr[3, 1], digits=2), '*horsepower) + (', 
                     round(leastSqr[4, 1], digits=2), '*weight) + (',           
                     round(leastSqr[5, 1], digits=2), '*acceleration) + (', 
                     round(leastSqr[1, 1], digits=2), ')')
## [1] "mpg = -0.01*displacement + (-0.04*horsepower) + (-0.01*weight) + (-0.02*acceleration) + (45.25)"

Reference:

https://github.com/wwells/CUNY_DATA_605/blob/master/Week05/WWells_Assign5.Rmd