echo on
%%%%% ch03_exa.m %%%%%
% Run script and study the effects of each line.

echo off
format compact;

format short;
disp(sprintf('**** Examples from Chapter 03 ****'));

disp(sprintf('**** Example 3.2&3.3 Normal Equations Method *'));
echo on
%%%% Linear lest square fitting by quadratic polynomial to points (ti,yi)
t=[-1.0 -0.5 0.0 0.5 1.0];	% five data points
y=[ 1.0  0.5 0.0 0.5 2.0]; 	% measured values 
A=[[ones(1,5)]' t' t'.^2]	% Vandermonde matrix for quadratic polynomial (5x3)
b=y';			% Right-hand side vector is built by measured values.

%%% Normal equation method A'Ax=A'b
Cp=A'*A		% cross-product matrix (nonsingular, symmetric, positive definite)(3x3)
bm=A'*b			% modified right-hand side
L=(chol(Cp))'		% Cholesky factorization gives upper triangular matrix L' (LL'x=bm).
y=L\bm;			% The solution of the lower triangular system.
x=L'\y			% The solution of the upper triangular system gives 
			% the polynomial coefficients that approximates (fits) given data.
         
         
disp(sprintf('**** Example 3.5 Householder QR Factorization *'));
echo on
%%%% Take matrix A and right-hand side vector b as in the previous example.
A		% Vandermonde matrix for quadratic polynomial (5x3)
bm=b		% Right-hand side vector is built by measured values.
R=A;		% Storage for R factor of A
I=eye(5);	% Identity matrix
H=eye(5);	% Storage for Householder transformation matrix
v=zeros(5,1);	% Temporary storage for vector v

%%%% annihilate the subdiagonal elements of the first column
a=R(:,1)%/max(R(:,1))	% Scale column R(:,1) with the largest element to prevent 
			% unnecessary overflow or underflow.
n=norm(a)		% calculate norm and
v=a;
if (a(1)<0) v(1)= v(1)-n; 	% take the same sign as R(1,1) to prevent cancellation
else v(1)= v(1)+n; end
v
H1=I-2*v*v'/(v'*v)	% The elementary elimination matrix H1
H=H1*H;		% The Householder transformation matrix H=H_n*H_n-1...H_2*H_1.
		% The elementary transformation matrices are not needed in general.
		% The Householder transformation can be implemented for each column of A separately
		% by the application of the vector v (complexity of O(3n))
R(:,1) = R(:,1)-((2*v'*R(:,1))/(v'*v))*v;	% on the first column of A
R(:,2) = R(:,2)-((2*v'*R(:,2))/(v'*v))*v;	% on the second column of A
R(:,3) = R(:,3)-((2*v'*R(:,3))/(v'*v))*v;	% on the third column of A
R
bm = bm -((2*v'*bm)/(v'*v))*v	% and on the right-hand side b

%%%% annihilate the subdiagonal elements of the second column
a=R(:,2)%/max(R(:,2))	% Scale column R(:,2) with the largest element to prevent 
			% unnecessary overflow or underflow.
a(1)=0;		% set all elements above diag. to 0 (annihilate only elements below diagonal).
n=norm(a)		% calculate norm and
v=a;
if (a(2)<0) v(2)= v(2)-n; 	% take the same sign as R(2,2) to prevent cancellation
else v(2)= v(2)+n; end
v
H2=I-2*v*v'/(v'*v)	% The elementary elimination matrix H2
H=H2*H;		% The Householder transformation matrix H=H_n*H_n-1...H_2*H_1.
		% The elementary transformation matrices are not needed in general.
		% The Householder transformation can be implemented for each column of A separately
		% by the application of the vector v (complexity of O(3n))
R(:,2) = R(:,2)-((2*v'*R(:,2))/(v'*v))*v;	% on the second column of A
R(:,3) = R(:,3)-((2*v'*R(:,3))/(v'*v))*v;	% on the third column of A
R
bm = bm -((2*v'*bm)/(v'*v))*v	% and on the right-hand side b

%%%% annihilate the subdiagonal elements of the last column
a=R(:,3)%/max(R(:,3))	% Scale column R(:,3) with the largest element to prevent 
			% unnecessary overflow or underflow.
a(1)=0; a(2)=0;	% set all elements above diag. to 0 (annihilate only elements below diagonal).
n=norm(a)		% calculate norm and
v=a;
if (a(3)<0) v(3)= v(3)-n; 	% take the same sign as R(2,2) to prevent cancellation
else v(3)= v(3)+n; end
v
H3=I-2*v*v'/(v'*v)	% The elementary elimination matrix H2
H=H3*H;		% The Householder transformation matrix H=H_n*H_n-1...H_2*H_1.
		% The elementary transformation matrices are not needed in general.
		% The Householder transformation can be implemented for each column of A separately
		% by the application of the vector v (complexity of O(3n))
R(:,3) = R(:,3)-((2*v'*R(:,3))/(v'*v))*v;	% on the third column of A
R
bm = bm -((2*v'*bm)/(v'*v))*v	% and on the right-hand side b
H			% The Householder transformation matrix is orthogonal, preserve the norm.
norm(H)
norm(bm)
norm(H*bm)
H*A		% Should be the same as R.

%%% Solution of the system Rx=bm 
bm=bm(1:3)		% Using only the first three components of bm
R=R(1:3,:)		% and upper triangular part of factorized A
x=R\bm		% gives already known solution for polynomial coefficients.
         
         
disp(sprintf('**** Example 3.7 Givens QR Factorization *'));
echo on
%%%% Take matrix A and right-hand side vector b as in the previous example.
A		% Vandermonde matrix for quadratic polynomial (5x3)
bm=b		% Right-hand side vector is built by measured values.
R=A;		% Storage for R factor of A
I=eye(5);	% Identity matrix
G=eye(5);	% Storage for Givens transformation matrix

%%%% annihilate the subdiagonal elements of the first column
%% start with bottom element (5,1)
Sr=sqrt(R(4,1)^2+R(5,1)^2); % Use (4,1)(also some other column element will work) and (5,1)
c=R(4,1)/Sr;		% appropriate rotations c 
s=R(5,1)/Sr;		% and s
Gi=I;
Gi(4,4)=c; Gi(5,5)=c;
Gi(5,4)=-s; Gi(4,5)=s;
G=Gi*G;			% Accumulate complete Givens transformation if needed.
R=Gi*R			% Element (5,1) is annihilated.
bm=Gi*bm		% Right hand side is modified.
%% proceed with the annihilation of the element (4,1)
Sr=sqrt(R(3,1)^2+R(4,1)^2); % Use (3,1) and (4,1)
c=R(3,1)/Sr;		% appropriate rotations c 
s=R(4,1)/Sr;		% and s
Gi=I;
Gi(3,3)=c; Gi(4,4)=c;
Gi(4,3)=-s; Gi(3,4)=s;
G=Gi*G;			% Accumulate complete Givens transformation if needed.
R=Gi*R			% Element (4,1) is annihilated.
bm=Gi*bm		% Right hand side is modified.
%%And so on....
echo off;
%% proceed with the annihilation of the element (3,1)
Sr=sqrt(R(2,1)^2+R(3,1)^2); % Use (2,1) and (3,1)
c=R(2,1)/Sr;		% appropriate rotations c 
s=R(3,1)/Sr;		% and s
Gi=I;
Gi(2,2)=c; Gi(3,3)=c;
Gi(3,2)=-s; Gi(2,3)=s;
G=Gi*G;			% Accumulate complete Givens transformation if needed.
R=Gi*R;			% Element (3,1) is annihilated.
bm=Gi*bm;		% Right hand side is modified.
%% proceed with the annihilation of the element (2,1)
Sr=sqrt(R(1,1)^2+R(2,1)^2); % Use (1,1) and (2,1)
c=R(1,1)/Sr;		% appropriate rotations c 
s=R(2,1)/Sr;		% and s
Gi=I;
Gi(1,1)=c; Gi(2,2)=c;
Gi(2,1)=-s; Gi(1,2)=s;
G=Gi*G;			% Accumulate complete Givens transformation if needed.
R=Gi*R;			% Element (2,1) is annihilated.
bm=Gi*bm;		% Right hand side is modified.
%% proceed with second column and element (5,2).
Sr=sqrt(R(4,2)^2+R(5,2)^2); % not to introduce nonzero elements take (4,2) and (5,2)
c=R(4,2)/Sr;		% appropriate rotations c 
s=R(5,2)/Sr;		% and s
Gi=I;
Gi(4,4)=c; Gi(5,5)=c;
Gi(5,4)=-s; Gi(4,5)=s;
G=Gi*G;			% Accumulate complete Givens transformation if needed.
R=Gi*R;			% Element (5,2) is annihilated.
bm=Gi*bm;		% Right hand side is modified.
%% proceed with the annihilation of the element (4,2).
Sr=sqrt(R(3,2)^2+R(4,2)^2); % take (3,2) and (4,2)
c=R(3,2)/Sr;		% appropriate rotations c 
s=R(4,2)/Sr;		% and s
Gi=I;
Gi(3,3)=c; Gi(4,4)=c;
Gi(4,3)=-s; Gi(3,4)=s;
G=Gi*G;			% Accumulate complete Givens transformation if needed.
R=Gi*R;			% Element (4,2) is annihilated.
bm=Gi*bm;		% Right hand side is modified.
%% proceed with the annihilation of the element (3,2).
Sr=sqrt(R(2,2)^2+R(3,2)^2); % take (2,2) and (3,2)
c=R(2,2)/Sr;		% appropriate rotations c 
s=R(3,2)/Sr;		% and s
Gi=I;
Gi(2,2)=c; Gi(3,3)=c;
Gi(3,2)=-s; Gi(2,3)=s;
G=Gi*G;			% Accumulate complete Givens transformation if needed.
R=Gi*R;			% Element (3,2) is annihilated.
bm=Gi*bm;		% Right hand side is modified.
%% proceed with the last column and element (5,3).
Sr=sqrt(R(4,3)^2+R(5,3)^2); % take (4,3) and (5,3)
c=R(4,3)/Sr;		% appropriate rotations c 
s=R(5,3)/Sr;		% and s
Gi=I;
Gi(4,4)=c; Gi(5,5)=c;
Gi(5,4)=-s; Gi(4,5)=s;
G=Gi*G;			% Accumulate complete Givens transformation if needed.
R=Gi*R;			% Element (5,3) is annihilated.
bm=Gi*bm;		% Right hand side is modified.

echo on;
%% proceed with the last column and element (4,3).
Sr=sqrt(R(3,3)^2+R(4,3)^2); % take (3,3) and (3,3)
c=R(3,3)/Sr;		% appropriate rotations c 
s=R(4,3)/Sr;		% and s
Gi=I;
Gi(3,3)=c; Gi(4,4)=c;
Gi(4,3)=-s; Gi(3,4)=s;
G=Gi*G			% Accumulate complete Givens transformation if needed.
R=Gi*R			% Element (4,3) is annihilated, final factorized A.
bm=Gi*bm		% Right hand side is modified.

norm(A)			% Norm of the original matrix is the same 
norm(R)			% as the norm of factorized matrix.
norm(b)			% Norm of the right-hand side vector is the same 
norm(bm)		% as the norm of modified vector.
Q=G'			% G*A=R and Q*G=I, because G is orthogonal.

%%% Solution of the system Rx=bm 
bm=bm(1:3)		% Using only the first three components of bm
R=R(1:3,:)		% and upper triangular part of factorized A
x=R\bm			% gives already known solution for polynomial coefficients.

disp(sprintf('**** Example 3.8 Modified Gram-Schmidt QR Factorization *'));
echo on
%%%% Take matrix A and right-hand side vector b as in the previous example.
A		% Vandermonde matrix for quadratic polynomial (5x3)
bm=b		% Right-hand side vector is built by measured values.
R=zeros(3);		% Storage for R factor of A
GS=A;		% Storage for Gram-Schmidt transformation matrix

%%%% Orthogonalize the first column against the subsequent columns
R(1,1)=norm(GS(:,1));	% Calculate the norm of the first column (R(1,1)) of A.
q1=GS(:,1)/R(1,1);		% Normalize the first column of A
R(1,2)=q1'*GS(:,2);	% Orthogonalizing the first column against the subsequent columns,
R(1,3)=q1'*GS(:,3);
GS=[q1, GS(:,2)-R(1,2)*q1, GS(:,3)-R(1,3)*q1] % gives the transformed matrix.

%%%% Orthogonalize the second column against the subsequent columns
R(2,2)=norm(GS(:,2));	% Calculate the norm of the second column (R(2,2)) of transformed A
q2=GS(:,2)/R(2,2);		% Normalize the second column of transformed A
R(2,3)=q2'*GS(:,3);	% Orthogonalizing the second column against remaining third column,
GS=[q1, q2, GS(:,3)-R(2,3)*q2] % gives the transformed matrix.

%%%% There is no more subsequent columns.
R(3,3)=norm(GS(:,3));	% Calculate the norm of the last column (R(3,3)) of transformed A
q3=GS(:,3)/R(3,3);		% Normalize the second column of transformed A
GS=[q1, q2, q3]; % gives the transformed matrix GS=Q 
Q=GS
R		% and R after Gram-Schmidt factorization.
Q*R	% So Q*R=A and the system Ax=b becomes QRx=b or Rx=bm
bm=Q'*b		% where bm is modified right-hand side,
x=R\bm	% with the solution as before.