Problem Set 1

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

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

(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.

Solving for \(\hat{x}\) in \(A^T A\hat{x} = A^T b\)

(x_hat <- inv(ATA)%*%ATb)
##      [,1]
## [1,]    1
## [2,]    4

What is the squared error of this solution?

Solving for \(e\) in \(e = b - A\hat{x}\)

# Error
(e <- b - A %*% x_hat)
##      [,1]
## [1,]   -1
## [2,]    3
## [3,]   -5
## [4,]    3
# Sum of squared error
sum(e^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 <- A%*%x_hat

#Zero vector to compare e
zero_v <- matrix(rep(0, length(p)), ncol = 1)
 
x2_hat <- inv(t(A) %*% A) %*% t(A) %*% p
(e2 <- p - A %*% x2_hat)
##               [,1]
## [1,]  8.881784e-16
## [2,]  0.000000e+00
## [3,]  0.000000e+00
## [4,] -3.552714e-15
(all.equal(zero_v, e2))
## [1] TRUE

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

e_actual <- matrix(c(-1,3,-5,3), ncol = 1)
all.equal(e_actual, b-p)
## [1] TRUE

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

If \(e\) is orthogonal to \(p\) and \(A\), their dot product will be zero.

# p orthogonal to e
all.equal(matrix(0), t(p) %*% e)
## [1] TRUE
# A orthogonal to e
all.equal(matrix(c(0,0)), t(A) %*% e)
## [1] TRUE

Problem Set 2

Using the auto-mpg data (linked from a repo - no need to change working directory) preform a least squares approximation to best fit the mpg in terms of the other data

mpg <- as.matrix(read.table(mpg_url))
A <- mpg[,-ncol(mpg)]
b <- mpg[,ncol(mpg)]

# Best fit solution 
(x_hat <- inv(t(A) %*% A) %*% t(A) %*% b)
##            [,1]
## V1 -0.030037938
## V2  0.157115685
## V3 -0.006217883
## V4  1.997320955
e <- b - A %*% x_hat

# Square of fitting error
(sum(t(e) %*% e))
## [1] 13101.43
#Plot of residuals
plot(e ~ jitter(A%*%x_hat),
     main = "Error from Least Squares Approximation",
     xlab = "A * x_hat",
     ylab = "Residuals")

abline(h = 0, lty = 3)