%% talked about the use of QR factorization for getting the least-squares
%% solution of an overdetermined linear system, as per notes on web.
% Stressed that condition number of A'*A is the square of that of A,
% while the condition number of the right-triangular factor in a qr
% factorization for A has the same condition number as A.

>> x = linspace(999,1001,11);
>> A = [x.^3 x.^2 x ones(size(x))];

>> size(A)
ans = 1 44
% I forgot to make x a column matrix!

>> x = linspace(999,1001,11)';
>> A = [x.^3 x.^2 x ones(size(x))];
>> cond(A)
ans =
  5.2747e+018
% get the economy size qr factorization:
>> [q,r] = qr(A,0);
>> r
r =
  1.0e+009 *
   -3.3166   -0.0033   -0.0000   -0.0000
         0   -0.0000   -0.0000   -0.0000
         0         0    0.0000    0.0000
         0         0         0    0.0000
>> invr = inv(r)
Warning: Matrix is close to singular or badly scaled.
         Results may be inaccurate. RCOND = 1.890167e-019.
invr =
  1.0e+009 *
   -0.0000    0.0000    0.0000   -0.0000
         0   -0.0000   -0.0000    0.0000
         0         0    0.0000   -0.0048
         0         0         0    1.5904

%% why is the inverse of an upper triangular matrix upper triangular???

>> cond(r)
ans = 5.2747e+018
>> cond(A)
ans = 5.2747e+018
>> cond(A'*A)
ans = 5.0580e+032
%% well, A is so badly conditioned, even the calculation of cond(A'*A) is
%% inaccurate
%  now choose a better way of writing cubic polynomials on the interval
%  [999 .. 1001], by using shifted powers:
>> m = mean(x)
m = 1000.0000
>> xm = x-m;
>> A = [xm.^3 xm.^2 xm ones(size(x))];
>> cond(A)
ans = 7.1036
>> cond(A'*A)
ans = 50.4606
>> [q,r] = qr(A,0);
>> cond(r)
ans =
    7.1036
%% ... i.e., exactly the same as cond(A)
%%% check its square, to compare with cond(A'*A):
>> ans^2
ans = 50.4606

%% ... exactly as expected.
>> diary off