Worksheet: Least squares approximation

In many applied scenarios, it is not practical (or even, perhaps, possible) to find an exact solution to a problem. Sometimes we may be working with imperfect data. Other times, we may be dealing with a problem that is overdetermined, ¹ such as a system of linear equations with more equations than variables. (Quite often, both of these issues may be present.)

🔗

An overdetermined system is quite likely to be inconsistent. That is, our problem requires finding a solution to a system

A x = b,

where no such solution exists! When no solution is possible, we can ask whether it is instead possible to find a best approximation.

🔗

What would a best approximation look like in this case? Let

U = col (A)

denote the column space of

A,

which we know is a subspace of

R^{n}

(assuming

A

is an

m \times n

matrix). The subspace

U

is precisely the set of all vectors

y

such that

A x = y

has a solution. Among these vectors, we would like to find the one that is closest to

b,

in the sense that

‖ y - b ‖

is as small as possible.

🔗

But we know exactly what this vector

y

should be: the orthogonal projection of

b

onto

U .

🔗

The presentation in this worksheet is based on the one given in the text by Nicholson (see Section 5.6). Further details may be found there, or, for the more statistically inclined, on Wikipedia ².

🔗

Given an inconsistent system

A x = b,

we have two problems to solve:

Find the vector $y = {proj}_{U} b,$ where $U = col (A)$
Find a vector $z$ such that $A z = y$

🔗

The vector

z

is then our approximate solution.

🔗

This can be done directly, of course:

Find a basis for $U$
Use the Gram-Schmidt Orthonormalization Algorithm to construct an orthogonal basis for $U$
Use this orthogonal basis to compute $y = {proj}_{U} b$
Solve the system $A z = y .$

🔗

But there is another way to proceed: we know that

b - y = b - A z \in U^{⊥},

so for any vector

A x \in U,

(A x) \cdot (b - A z) = 0 .

Therefore,

\begin{aligned} (A x) \cdot (b - A z) & = 0 \\ (A x)^{T} (b - A z) & = 0 \\ x^{T} A^{T} (b - A z) & = 0 \\ x^{T} (A^{T} b - A^{T} A z) & = 0, \end{aligned}

🔗

for any vector

x \in R^{n} .

🔗

Therefore,

A^{T} b - A^{T} A z = 0,

or

\begin{matrix} (3.5.1) & A^{T} A z = A^{T} b . \end{matrix}

🔗

Solving this system, called the normal equations for

z,

will yield our approximate solution.

🔗

To begin, let’s compare the two methods discussed above for finding an approximate solution. Consider the system of equations

A x = b,

where

A = [\begin{matrix} 3 & - 1 & 0 & 5 \\ - 2 & 7 & - 3 & 0 \\ 4 & - 1 & 2 & 3 \\ 0 & 3 & 9 & - 1 \\ 7 & - 2 & 4 & - 5 \\ 1 & 0 & 3 & - 8 \end{matrix}] and b = [\begin{matrix} 4 \\ 0 \\ 1 \\ 2 \\ - 5 \\ - 1 \end{matrix}] .

🔗

1.

🔗

Confirm that the system has no solution.


    
        
xxxxxxxxxx
 
1
from sympy import Matrix,init_printing
2
init_printing()

    
    
    
    
        
            
                Language:
                
            
        
    
    




    
    
        
        Messages


    
        
xxxxxxxxxx
 
1
In Jupyter, double-click to edit, and change this to a markdown cell to explain your results.

    
    
    
    
        
            
                Language:
                
            
        
    
    




    
    
        
        Messages

🔗

2.

🔗

Find an orthogonal basis for the column space of

A .


    
        
xxxxxxxxxx
 
1
from sympy import GramSchmidt

    
    
    
    
        
            
                Language:
                
            
        
    
    




    
    
        
        Messages

🔗

3.

🔗

Compute the projection

y

of

b

onto the column space of

A .

🔗

4.

🔗

Solve the system

A z = y

for

z .

🔗

5.

🔗

Solve the normal equations (3.5.1) for

z .

🔗

Next, we want to consider a problem found in many introductory science labs: finding a line of best fit. The situation is as follows: in some experiment, data points

(x_{1}, y_{1}), (x_{2}, y_{2}), \dots, (x_{n}, y_{n})

have been found. We would like to find a function

y = f (x) = a x + b

such that for each

i = 1, 2, \dots, n,

the value of

f (x_{i})

is as close as possible to

y_{i} .

🔗

Note that we have only two parameters available to tune:

a

and

b .

We assume that some reasoning or experimental evidence has led us to conclude that a linear fit is reasonable. The challenge here is to make precise what we mean by “as close as possible”. We have

n

differences (sometimes called residuals)

r_{i} = y_{i} - f (x_{i})

that we want to make small, by adjusting

a

and

b .

But making one of the

r_{i}

smaller might make another one larger!

🔗

A measure of the overall error in the fit of the line is given by the sum of squares

S = r_{1}^{2} + r_{2}^{2} + \dots + r_{n}^{2},

🔗

and this is the quantity that we want to minimize. (Hence the name, “least squares”.)

🔗

Let

v = [\begin{matrix} a \\ b \end{matrix}],

and note that

f (x) = a + b x = [\begin{matrix} 1 & x \end{matrix}] v .

Set

y = [\begin{matrix} y_{1} \\ y_{2} \\ ⋮ \\ y_{n} \end{matrix}]

and

r = [\begin{matrix} r_{1} \\ r_{2} \\ ⋮ \\ r_{n} \end{matrix}] .

Then

r = y - A v,

🔗

where

A = [\begin{matrix} 1 & x_{1} \\ 1 & x_{2} \\ ⋮ & ⋮ \\ 1 & x_{n} \end{matrix}] .

(Note that we are using

y

to denote a different sort of vector than on the previous page.)

🔗

We can safely assume that an exact solution

A v = y

is impossible, so we search for an approximate one, with

r

as small as possible. (Note that the magnitude of

r

satisfies

‖ r ‖^{2} = S .

) But a solution

z

that makes

y - A z

as small as possible is exactly the sort of approximate solution that we just learned to calculate! Solving the normal equations for

z = [\begin{matrix} a \\ b \end{matrix}],

we find that

z = (A^{T} A)^{- 1} (A^{T} y) .

🔗

6.

🔗

Find the equation of the best fit line for the following set of data points:

(1, 2.03), (2, 2.37), (3, 2.91), (4, 3.58), (5, 4.11), (6, 4.55), (7, 4.93), (8, 5.44), (9, 6.18) .

🔗

7.

🔗

Suppose we were instead trying to find the best quadratic fit for a dataset. What would our parameters be? What would the matrix

A

look like? Illustrate with an example of your own.

Linear Algebra: A second course, featuring proofs and Python

Worksheet 3.5 Worksheet: Least squares approximation
A4 US

1.

2.

3.

4.

5.

6.

7.

Worksheet 3.5 Worksheet: Least squares approximation A4US

1.

2.

3.

4.

5.

6.

7.

Worksheet 3.5 Worksheet: Least squares approximation
A4 US