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);

>

>