%%  Example of  QR decomposition.   Application of QR to Least Squares.


clear, close all

% Create a random matrix

A=randn(4,3)

q1=A(:,1);
R11 = norm(A(:,1));
q1 = q1/R11;

R12 = q1'*A(:,2);
q2=A(:,2) - R12*q1;
R22 = norm(q2);
q2=q2/R22;

R13=q1'*A(:,3);
R23=q2'*A(:,3);
q3=A(:,3) - R13*q1-R23*q2;
R33=norm(q3);
q3 = q3/R33;

%% The GS yields the orthonormal vectors:

q1
q2
q3

%% Form the orthogonal matrix Q:

Q=[q1,q2,q3]

%% Form the R matrix

R=[R11, R12, R13; 0, R22, R23; 0, 0, R33]
fprintf('check A=QR by finding A-QR \n')
A-Q*R

%% A Small Application of QR to Least Squares.
%% We will solve for z, the solution of
%%  A z = Pb,
%%   where Pb is the projection of b onto col(A).  We create b as a
%% product of q1.*q1.*q1 (i.e. a direct product of q1 with itself. This will guarantee that the vector so produced will not be entirely in col(A).

b = q1.^3

%% We now use the A=QR
%%  to solve a problem Ax = b in the least squares sense:
%% The projection Pb = Q Q^T of b onto col(A):

Q*Q'*b

%% The least squares solution  R zz = Q^Tb. It is zz = R^{-1} Q^T b,
%%  but since R is upper triangular

R

%% and Q^T b

qtb=Q'*b;

%% we can quickly backsubstitute:

z3=qtb(3,1)/R33;
z2 =(-R23*z3+qtb(2,1))/R22;
z1 =(-R12*z2-R13*z3+qtb(1,1))/R11;
zz=[z1;z2;z3]

%%  Let's check: Azz-QQ^T b

A*zz-Q*Q'*b