```% MATLAB Demo (Least squares)
% File: lsdemo
%
% Find the least-squares parabola for the 5 data points:
%
%            t  |    y
%          ==============
%            0  |   -3
%            1  |    7.5
%            2  |    2.5
%            3  |   14.5
%            4  |   23.5
%
% Formulate problem in terms of finding the least-squares solution to Ax = b.
% Verify the results, and plot a graph of the data vs. the parabola.
% Compute the projection matrix for the column space of A.
%
%  Outline:
%  The matrix A has 5 rows; row k is [ 1   t_{k}   t_{k}^2 ].
%  The vector b has 5 components; component k  is y_{k}.
%  Use  ref([A b]) to show that Ax = b is NOT solvable.
%  Use 18.06 MATLAB command  lsq(A, b) to find the least-squares solution.
%  Verify: p = A*xbar, p+e = b, and e is orthogonal to p.
%  Give a plot of data points vs. the least-squares parabola.
%  Compute the projection matrix, or use 18.06 command  projmat(A)  .

>> diary lsdemo
>> b = [-3; 7.5; 2.5; 14.5; 23.5];
>> A = [1 1 1 1 1; 0 1 2 3 4; 0 1 4 9 16];
>> A = A'
A =
1     0     0
1     1     1
1     2     4
1     3     9
1     4    16

>> Z = [A b]
Z =
1.0000         0         0   -3.0000
1.0000    1.0000    1.0000    7.5000
1.0000    2.0000    4.0000    2.5000
1.0000    3.0000    9.0000   14.5000
1.0000    4.0000   16.0000   23.5000

>> Z = ref(Z)
Z =
1     0     0     0
0     1     0     0
0     0     1     0
0     0     0     1
0     0     0     0
%
%%% This reduced row echelon form [R d] shows that d cannot
%%% be expressed as a linear combination of the columns of R.
%%% Rx=d (and Ax=b) is not solvable.
%
>> help lsq
LSQ	Least squares.
[xbar,p,e] = LSQ(A,b) finds a least squares
solution to the overdetermined system A*x ~= b.
xbar = the solution to the normal equations.
p = projection of b onto the column space.
e = b - p.

>> [xbar, p, e] = lsq(A, b)
xbar =
-1
2
1
p =
-1
2
7
14
23
e =
-2.0000
5.5000
-4.5000
0.5000
0.5000

%%% Check: A*xbar = p.
>> A * xbar
ans =
-1
2
7
14
23

%%% Check: p+e = b.
>> p+e
ans =
-3.0000
7.5000
2.5000
14.5000
23.5000

%%% Check: e is orthogonal to p.
>> p'*e
ans =
0
%
%%% Plot the data points and parabola over the interval -1 <= t <= 5.
%
>> t = linspace(-1, 5, 21);
>> t = t';
>> xbar
xbar =
-1
2
1
%
%%% Compute y(t) for each value of t in the interval [-1, 5].
%
>> y = xbar(1) + xbar(2)*t + xbar(3)*t.^2;
>> plot(t,y,'r-');
>> grid on; hold on
%
%%% Plot the 5 data points, and label the axes.
%
>> plot([0; 1; 2; 3; 4], b, 'g+');
>> xlabel('t');
>> ylabel('y(t)');
>> title('Least-squares parabola: y(t)');
>> print -dps parabola.ps

%
%%% Compute the projection matrix, P = A*inv(A'*A)*A'.
%%% Recall that A has 5 rows and 3 columns.
%%% A must have linearly independent columns in order for A'*A
%%% to be invertible.
%
>> format rat
>> P = A*inv(A'*A)*A'
P =
31/35         9/35        -3/35        -1/7          3/35
9/35        13/35        12/35         6/35        -1/7
-3/35        12/35        17/35        12/35        -3/35
-1/7          6/35        12/35        13/35         9/35
3/35        -1/7         -3/35         9/35        31/35

%%% Check: the columns of P span a 3-dimensional subspace of R^{5}.
%%% In particular, this subspace is the column space of the matrix A.
>> rank(P)
ans =
3

>> help projmat
PROJMAT Projection matrix onto the column space.
P = projmat(A) is the symmetric, square matrix that
projects any vector onto the column space of A.

>> projmat(A)
ans =
31/35         9/35        -3/35        -1/7          3/35
9/35        13/35        12/35         6/35        -1/7
-3/35        12/35        17/35        12/35        -3/35
-1/7          6/35        12/35        13/35         9/35
3/35        -1/7         -3/35         9/35        31/35

%%% Check: The projection of b onto the column space of A is p.
>> projmat(A)*b
ans =
-1
2
7
14
23
>> diary off
```