__ ____ ____Least Squares - The Two Dimensional Case__

Solving an over-determined linear system A x = b is generally impossible, but

finding an x0 so that A x0 is as close as possible to b leads to the concept of a

least squares solution. One can see that such an x0 must satisfy the

normal equations

A_tr A x0 = A_tr b

where A_tr is the transpose of the matrix A. Generically this system will

have a unique solution. We use this worksheet to study this in the case of

fitting a set of points to a line.

**Problem: **
Find the best fitting straight line y = cx + d through the following points:

(1,2), (1,3),(2,4), (3,5)

`> `
**pointlist := [[1,2],[1,3],[2,4],[3,5]];**

**Solution:**
The over-determined system here is

1 c + d = 2

1 c + d = 3

2 c + d = 4

3 c + d = 5

In matrix form: A x = b:

`> `
**with(linalg):**

Warning, the protected names norm and trace have been redefined and unprotected

`> `
**A := array([[1,1],[1,1],[2,1],[3,1]]);**

`> `
**x :=vector([c,d]);**

`> `
**b := vector([2,3,4,5]);**

This step just confirms we've chosen A and b to have the right over-determined system.

`> `
**evalm(A &* x = b );**

The normal equations N x = m :

`> `
**N := evalm(transpose(A) &* A);**

`> `
**m := evalm(transpose(A) &* b);**

`> `
**x0 := evalm(inverse(N) &* m);**

`> `
**line := x -> x0[1] * x + x0[2];**

`> `
**plot({pointlist,line},
title = `Least Squares Solution Vs. Points Connected`);**

Here's a singular example where the normal equations don't have a unique solution.

`> `
**pointlist := [[1,1],[1,2],[1,3],[1,4]];**

`> `
**A:= array([[1,1],[1,1],[1,1],[1,1]]);**

`> `
**b := vector([1,2,3,4]);**

The normal equations N x = m :

Since A had rank 1, N also has rank 1 and is not invertible.

`> `
**N := evalm(transpose(A) &* A);**

Still the equations N x = m have solutions - the normal equations always do !

`> `
**m := evalm(transpose(A) &* b);**

It is tempting to investigate what happens when the problem is nearly singular.

We change the first point to (1.000001,1.00001).

`> `
**A:= array([[1.000001,1],[1,1],[1,1],[1,1]]);**

`> `
**b := vector([1.000001,2,3,4]);**

`> `
**N := evalm(transpose(A) &* A);**

`> `
**m := evalm(transpose(A) &* b);**

`> `
**x0 := evalm(inverse(N) &* m);**

Error, (in inverse) singular matrix

`> `
**Digits;**

`> `
**Digits := 40;**

`> `
**x0 := evalm(inverse(N) &* m);**

`> `

`> `