clear
echo on
%In this example we demonstrate how
%QR factorization can be used to solve
%a linear regression problem.  First
%let's set up the problem.  Suppose we
%want to measure the acceleration due to
%gravity by measuring the position of a
%ball as a function of the time after it
%has been dropped.  We know that the
%position should be given by:
%
%   h = h0 + u0*t -0.5*g*t^2
%
%where h0 is the initial position and u0
%is the initial velocity.
%Thus if we have the data points (t(i),h(i))
%we can fit a parabola to it and get g.
pause
%We take the data to be given by:
t=[0,.1,.2,.3,.4,.5,.6]'
%and
h=[0.5,14.6,21.0,15.5,2.0,-22.7,-56]'
pause
%We can plot this data:
plot(t,h,'o')
xlabel('time')
ylabel('height')
pause
%The modeling functions are 1, t, and t^2,
%thus the matrix A is given by:
n=length(t);
a=[ones(n,1),t,t.^2]
pause
%We can solve this using the normal equations:
ata=a'*a
atb=a'*h
pause
%and then the solution:
x=ata\atb
%with the gravitational acceleration of:
g=-2*x(3)
pause
%We can compare this to the data:
hold on
plot(t,a*x)
pause
%Now let's do the same thing via QR factorization.
%First we do the QR factorization:
[q,r]=qr(a);
%Let's look at r first:
r
%Note that only upper triangular elements are non-
%zero.  The diagonal elements of this matrix are
%not necessarily ordered by size, although 'qr' can be set
%using pivoting to be so.
pause
%Next we look at the matrix q:
q
%This is a square matrix with the interesting property
%that q'*q=eye, that is q is orthogonal.
q'*q
pause
%We can reduce the problem by multiplying both A and
%the rhs h by q':
qta=q'*a
%and
qtb=q'*h
pause
%Note that the last four elements of the new rhs
%provide the norm of the residual.
%OK, let's solve this:
x=r\qtb
%which gives the same result.
echo off