% CP3p4 on p.106 of Heath
% Example of a linear least squares problem that is sensitive.
% Solve the normal equations using Cholesky factorization.
% K. Ming Leung, 04/07/03

clear all; format long g;
AA = [ 0.16 0.10
       0.17 0.11
       2.02 1.29];

% 2 slightly different right-hand side vectors
Bvec1 = [ 0.26; 0.28; 3.31]; 
Bvec2 = [0.27; 0.25; 3.33];
 
n = min(size(AA));
AT = AA.';
ATA = AT*AA;             
ATb1 = AT*Bvec1;
ATb2 = AT*Bvec2;
% Solve the system of normal equations using Cholesky factorization
LL = Cholesky(ATA); 
Vec1 = forwardSubstitution(LL,ATb1,n);
Xvec1 = backSubstitution(transpose(LL),Vec1,n),
Vec2 = forwardSubstitution(LL,ATb2,n);
Xvec2 = backSubstitution(transpose(LL),Vec2,n),
% compute using Matllab's backslash operator for comparison
Xmatlab1 = AA\Bvec1,
Xmatlab2 = AA\Bvec2,

conditionA = cond(AA),

Rvec1 = Bvec1-AA*Xvec1, normr1 = sqrt(Rvec1.'*Rvec1), 
Rvec2 = Bvec2-AA*Xvec2, normr2 = sqrt(Rvec2.'*Rvec2), 
normAx1 = norm(AA*Xvec1), normb1 = norm(Bvec1), cosTheta1 = normAx1/normb1,
normAx2 = norm(AA*Xvec2), normb2 = norm(Bvec2), cosTheta2 = normAx2/normb2,

% The most striking result is that even though b1 and b2 
% differ by only about 4%, the solutions differ by several 100%.
% Note that the solutions obtained by solving the set of normal equations
% agree very well those obtained from Matlab (which uses QR factorization).
% So the problem is not due to the use of normal equations.
%
% The problem is ill-conditioned since
% matrix A has a condition number of about 1098.
% Note that the residual vectors are small in both cases.
% Thus the angle between b and y=A*x, theta, is small, and
% there is no condition squaring effect (see p.115 of textbook).
% According to the equation at the top of p.115,
% the relative error in the solution is given by
% cond(A)/cos(theta) times the relative error in b,
% and so is about 1098/1*0.04, which is about 4.
% This is in agreement with our numerical finding.
