\[A = \left[ \begin{array}{cc}1 & 0 \\1 & 1 \\1 & 3 \\1 & 4 \end{array} \right]\ ; \ \ \ b = \left[ \begin{array}{c}0 \\8 \\8 \\20 \end{array} \right]\]
A = matrix(c(1, 1, 1, 1, 0, 1, 3, 4), ncol = 2, byrow = FALSE)
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\[A^TA = \left[ \begin{array}{cc}4 & 8 \\8 & 26 \end{array} \right]\ ; \ \ \ A^Tb = \left[ \begin{array}{c}36 \\112 \end{array} \right]\]
(x_hat <- solve(AtA, Atb)) [,1]
[1,] 1
[2,] 4\[\hat{x} = \left[ \begin{array}{c}1 \\4 \end{array} \right]\]
(e <- b - A %*% x_hat) [,1]
[1,] -1
[2,] 3
[3,] -5
[4,] 3\[e = \left[ \begin{array}{c}-1 \\3 \\-5 \\3 \end{array} \right]\ ; \ \ \ ||e||^2 = 44\]
p <- matrix(c(1, 5, 13, 17), ncol = 1)
Atp <- t(A) %*% p
solve(AtA, Atp) [,1]
[1,] 1
[2,] 4\[\hat{x_p} = \left[ \begin{array}{c}1 \\4 \end{array} \right]\]
identical(e, b - p)[1] TRUEsum(e * p)[1] 0for (i in 1:ncol(A)) {print(sum(e * A[, i]))}[1] 0
[1] 0The dot products of \(e\) with \(p\) and with the columns of \(A\) yield zero; this shows that the vectors are orthogonal.
The function below returns the inverse of a full-rank square matrix using co-factors.
# read in data and store in A & b
car_data <- read.table('https://raw.githubusercontent.com/dsmilo/DATA605/master/data/auto-mpg.data',
col.names = c("x1", "x2", "x3", "x4", "mpg"))
A <- as.matrix(car_data[, 1:4])
b <- as.matrix(car_data[, 5])
# find AtA & Atb
AtA <- t(A) %*% A
Atb <- t(A) %*% b
# find least squares solution
(x_hat <- solve(AtA, Atb)) [,1]
x1 -0.030037938
x2 0.157115685
x3 -0.006217883
x4 1.997320955# find the error
e <- b - A %*% x_hat
e2 <- sum(e^2)The solution vector returned by this script is \[\hat{x} = \left[ \begin{array}{c}-0.030 \\0.157 \\-0.006 \\1.997 \end{array} \right]\]
That is to say that the least-squares estimate for \(mpg\) is \[\widehat{mpg} = -0.030 x_1 + 0.157 x_2 -0.006 x_3 + 1.997 x_4\]
The fitting error of this estimate, \(||e||^2\), is 13101.43.