% 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