Least Squares Problem

Least Squares Problem #

Implementing Least Squares problem solver, along with some examples of it running.

.rmd knitted with execution results available for download here

Define a function to solve the least squares problem#

Function takes in a list of (x,y) points.

Creates a matrix ‘A’, with the first column being all 1’s, and the second column being all the ‘x’ values

Creates a vector ‘b’, with it’s values being all of the ‘y’ values

Then solves the equation: x~ = inv(t(A)A)t(A)b

… To get x~, the ‘a’ and ‘b’ values of the y = a + bx line of best fit

1
least_squares_linear = function(points){
2
  # create 'A'
3
  A = matrix(c(rep(1, length(points$x)), points$x), nr=length(points$x), nc=2)
4
  # create 'b'
5
  b = points$y
6
  # calculate 't(A)A'
7
  AtA = t(A) %*% A
8
  AtA
9
  # find the inverse of 'AtA'
10
  invAtA = solve(AtA)
11
  # use the previous values to solve for 'x~'
12
  solution_linear = invAtA %*% t(A) %*% b
13

14
  return(solution_linear)
15
}

The only difference for the quadratic version is the additional column in ‘A’ containing all the ‘x’ values, squared

Realistically, this could be worked into a general function that works for polynomials to the ’n’th degree, but I can’t be bothered.

1
least_squares_quad = function(points){
2
  # create 'A'
3
  A = matrix(c(rep(1, length(points$x)), points$x, points$x ^ 2), nr=length(points$x), nc=3)
4
  # create 'b'
5
  b = points$y
6
  # calculate 't(A)A'
7
  AtA = t(A) %*% A
8
  AtA
9
  # find the inverse of 'AtA'
10
  invAtA = solve(AtA)
11
  # use the previous values to solve for 'x~'
12
  solution_quad = invAtA %*% t(A) %*% b
13

14
  return(solution_quad)
15
}

Problem 1:#

Find the least squares approximating line for the points: * (1,0) * (2,2) * (4,5) and compute the corresponding least squares error.

1
# least squares error is the difference between the actual point, and the predicted point, squared, summed, and then square rooted
2
error = function(a, b, points){
3
  e_i = points$y - (a + b*points$x)
4
  return(sqrt(sum(e_i^2)))
5
}
6

7
points = list()
8
points$x = c(1,2,4)
9
points$y = c(0,2,5)
10

11
sol = least_squares_linear(points)
12
# least squares solution is:
13
sol

1
##           [,1]
2
## [1,] -1.500000
3
## [2,]  1.642857

1
plot(x=points$x, y=points$y)
2
lines(x=points$x, y=sol[1] + sol[2] * points$x)

```r e = error(sol[1], sol[2], points) # least squares error is: e ``` ``` ## [1] 0.2672612 ```

Problem 2:#

Find the least squares approximating parabola for the points

(1,1)
(2,-2)
(3,3)
(4,4)

Since parabola’s are of degree 2, use the quadratic least squares function

1
points = list()
2
points$x = c(1,2,3,4)
3
points$y = c(1,-2,3,4)
4

5
sol = least_squares_quad(points)
6
# solution is:
7
sol

1
##      [,1]
2
## [1,]  3.0
3
## [2,] -3.6
4
## [3,]  1.0

1
plot(x=points$x, y=points$y)

1
curve(sol[1] + sol[2]*x + sol[3]*x^2, from=min(points$x), to=max(points$x))

Problem 3:#

Construct and solve the normal equations to find a least squares solution of Ax = b, where

1
A = matrix(c(c(3,1,1), c(1,2,1)), nc=2, nr=3)
2
A

1
##      [,1] [,2]
2
## [1,]    3    1
3
## [2,]    1    2
4
## [3,]    1    1

1
b = c(1,1,1)
2
b

1
## [1] 1 1 1

This is just solving using the linear least squares again, but we have ‘A’ and ‘b’ instead of points

So, just use the linear least squares function code to solve for ‘x’

1
# calculate 't(A)A'
2
AtA = t(A) %*% A
3
# find the inverse of 'AtA'
4
invAtA = solve(AtA)
5
# use the previous values to solve for 'x~'
6
sol = invAtA %*% t(A) %*% b
7

8
sol

1
##           [,1]
2
## [1,] 0.2000000
3
## [2,] 0.4666667

Problem 4:#

Consider

1
A = matrix(c(c(1,1,0,0),c(1,0,-1,-1),c(0,1,1,1),c(0,1,1,0)), nc=4, nr=4)
2
A

1
##      [,1] [,2] [,3] [,4]
2
## [1,]    1    1    0    0
3
## [2,]    1    0    1    1
4
## [3,]    0   -1    1    1
5
## [4,]    0   -1    1    0

1
b = c(1,-3,2,4)
2
b

1
## [1]  1 -3  2  4

Show that the least squares solution of Ax = b is not unique. Solve the normal equations to find all least squares solutions.

[Hint: The RREF might be helpful.]

1
if (!require(pracma)) {
2
  install.packages("pracma")
3
}

1
## Loading required package: pracma

1
library(pracma)
2

3
rref(A)

1
##      [,1] [,2] [,3] [,4]
2
## [1,]    1    0    1    0
3
## [2,]    0    1   -1    0
4
## [3,]    0    0    0    1
5
## [4,]    0    0    0    0

Theorem 14 in the course notes states that the least squares solution is unique if ‘A’ has linearly independent columns, but this is not the case, since the result of RREF(A) has a zero row.

Therefore the solution to the least squares problem is not unique.

1
# solving for a general solution
2
# find t(A)A
3
AtA = t(A) %*% A
4
AtA

1
##      [,1] [,2] [,3] [,4]
2
## [1,]    2    1    1    1
3
## [2,]    1    3   -2   -1
4
## [3,]    1   -2    3    2
5
## [4,]    1   -1    2    2

1
# find t(A)b
2
Atb = t(A) %*% b
3
Atb

1
##      [,1]
2
## [1,]   -2
3
## [2,]   -5
4
## [3,]    3
5
## [4,]   -1

1
# form augmented matrix, and row reduce
2

3
aug = matrix(c(AtA, Atb), nr=4, nc=5)
4
rref(aug)

1
##      [,1] [,2] [,3] [,4] [,5]
2
## [1,]    1    0    1    0    3
3
## [2,]    0    1   -1    0   -4
4
## [3,]    0    0    0    1   -4
5
## [4,]    0    0    0    0    0

1
# formulate the general solution using solved matrix
2
# shuffling the equations around and solving for each 'xi'...
3
# the free variable is 'x3', meaning it can be anything, for now it's just set to 0
4
x3 = 0
5
x1 = 3 - x3
6
x2 = -4 + x3
7
x4 = -4
8

9
# so the solution is:
10
x_tilde = function(x3){
11
  return(c(3,-4,0,-4) + x3 * c(-1, 1, 1, 0))
12
}
13
sol1 = x_tilde(0)
14
sol1

1
## [1]  3 -4  0 -4

1
sol2 = x_tilde(1)
2
sol2

1
## [1]  2 -3  1 -4