%% Week 5
%% Householder reflection, QR decomposition. Givens rotation and QR decomposition with Givens rotations. Solution of over-determined systems
%
%% QR decomposition
%
% In this topic, the matrices will be matrices with full column rank. This means that if the matrix has m rows and n columns then its rank is r\left(A\right)=n (the column vectors are linearly independent).
% Our goal is to decompose a matrix A\in \mathbb{R}^{m\times n}, r(A)=n in the form A=QR, where Q is an mxm orthogonal matrix (its inverse is its transpose) and R is an mxn upper triangular matrix.
% Later we will see two techniques to perform this decomposition but first let us see how Matlab computes Q and R.
% Example 1. (QR decomposition with Matlab)

A=[1 2 3
   4 5 6
   7 8 9
   10 11 15]

[Q,R]=qr(A) % The command is simply qr, the output parameters give the matrices Q and R.

% We can see that R is an upper triangular matrix and Q is ortohogonal
% because

Q'*Q % is the identity matrix

% and the product QR gives back the original matrix A.

Q*R

% The main idea in the construction of the QR decomposition is as follows. We multiply the matrix A by orthogonal matrices, say S_1 ,S_2 ,\ldotp \ldotp \ldotp ,S_k, from left until we obtain an upper triangular matrix from A. Formally:
% {S_{k\;} \ldotp \ldotp \ldotp \;S}_2 S_1 A=R\;\left(R\;\mathrm{is}\;\mathrm{already}\;\mathrm{an}\;\mathrm{upper}\;\mathrm{triangular}\;\mathrm{matrix}\right). Then we can write
% A=S_1^T S_2^T \ldotp \ldotp \ldotp S_k^T R=\mathrm{QR}, where the matrix S_1^T S_2^T \ldotp \ldotp \ldotp S_k^T is an orthogonal matrix (product of orhogonal matrices) and can be denoted by Q.
% Example 2. (QR decomposition with a sequence of orthogonal matrices)

A=[1 2 3
   4 5 6
   7 8 9
   10 11 15]

% Let us define the following matrices. 

S1=[-0.077615052570633  -0.310460210282533  -0.543305367994433  -0.776150525706333
  -0.310460210282533   0.910556611158364  -0.156525930472862  -0.223608472104089
  -0.543305367994433  -0.156525930472862   0.726079621672491  -0.391314826182156
  -0.776150525706333  -0.223608472104089  -0.391314826182156   0.440978819739778]

S2=[   1.000000000000000                   0                   0                   0
                   0  -0.211234790319991   0.350582131484409   0.912399053288811
                   0   0.350582131484409   0.898526832371054  -0.264086539969549
                   0   0.912399053288811  -0.264086539969549   0.312707957948936]

S3=[   1.000000000000000                   0                   0                   0
                   0   1.000000000000000                   0                   0
                   0                   0  -0.997688645856591   0.067951202556251
                   0                   0   0.067951202556251   0.997688645856592]

% These matrices are orhogonal. Let us check this!

S1*S1'
S2*S2'
S3*S3' % These are identity matrices indeed.

% Construct the product:

R=S3*S2*S1*A % This is an upper triangular matrix that can be denoted by R.

% Thus, matrix Q can be chosen to be

Q=S1'*S2'*S3'

% We have obtained the same Q and R matrices as in the previous example.

% Of course, the question is how to construct the matrices S_1 ,S_2 ,\ldotp \ldotp \ldotp that can transform matrix A into an upper triangular form. There are two types of matrices that do this job: Householder reflections and Givens rotations.
%% Householder reflection
%
% It is easy to see that the matrix H=I-2\frac{vv^T}{v^Tv}  (I is the identity matrix) reflects the vector x through the plane that contains the origin and is orthogonal to the vector v (v is the normal vector of the reflection plane). This fact is based on the form
% {x}^\prime={x}-\frac{2{x}^T{x}}{{x}^T{x}}v=x-\frac{2vv^Tx}{v^Tv}=(I-\frac{2vv^T}{v^Tv})x
% of the reflected vector. This figure demonstrates the reflection in 2D.
%
% Example 3. (reflection through a plane)

% What is the reflection of the vector x=[1;2;4] through the plane with
% normal vector v=[3;4;5]?

x=[1;2;4]
v=[3;4;5]

H=eye(3)-2*v*v'/(v'*v)

x_reflected=H*x % This will be the reflected image.

% Just to see the vectors and the plane in 3D:

figure
plot3([0,x(1)],[0,x(2)],[0,x(3)])
axis equal
hold on
plot3([0,x_reflected(1)],[0,x_reflected(2)],[0,x_reflected(3)])
hold on
[X,Y]=meshgrid([-4:0.1:2]);
mesh(X,Y,(-v(1)*X-v(2)*Y)/v(3))
hold off
snapnow

%
% Two important properties of the reflection matrix H are formulated in Theorem 34.
% Thm. 34. H is a symmetric and orthogonal (we need this property to the QR decomposition) matrix.
% Proof. The symmetry is trivial. Moreover
% \left(I-\frac{2vv^T}{v^Tv}\right)\left(I-\frac{2vv^T}{v^Tv}\right)=I-4\frac{vv^T}{v^Tv}+4\frac{vv^T}{v^Tv}\frac{vv^T}{v^Tv}=I.
%
% Because we would like to use reflection matrices to construct the QR decomposition, it would be nice to choose them (this practically means the choice of the normal vector v) such a way that they zero out as many elements of x as possible. This is achievable. We can zero out all elements of a vector except the first one by choosing the normal vector appropriately. We can reflect the vector x to the first coordinate axis.
%
%
%
%
% Example 4. (Householder reflection)

% Let us give the Householder reflection matrix to the vector x=[-3;4;3;2]
% and perform the transformation.

x=[-3;4;3;2]

% We determine the normal vector of the reflection plane

v=x;
v(1)=v(1)-norm(x) % This will be the normal vector. Here we used subtraction because 
% the first element of x is negative (to avoid catastrophic cancellation)

% We constuct the Householder reflection

H=eye(4)-2*v*v'/(v'*v)

% The result of the transformation

H*x % Only the first element is nonzero. Its value is the norm of the vector x 
% (reflection does not change the 2-norm of the vectors).

% With Householder reflection we can construct the QR decomposition of a matrix. This procedure is used also to prove the existence of the QR decomposition of full rank matrices.
%
% Example 5. (The QR decomposition of matrix A with Householder reflections from Example 1.)

A=[1 2 3
   4 5 6
   7 8 9
   10 11 15]

% The first Householder matrix to the first column.

v=A(:,1);
v(1)=v(1)+norm(v);
H1=eye(4)-2*v*v'/(v'*v)

% this zeros out the first column of A below the diagonal

A1=H1*A

% The second Householder matrix to the second column

v=A1(2:4,2);
v(1)=v(1)+norm(v);
H2=blkdiag(1,eye(3)-2*v*v'/(v'*v)) % this is a blockdiagonal matrix. The (1,1)
% element is 1 and the (2:4,2:4) block is the Householder reflection.

% this zeros out the second column below the main diagonal

A2=H2*H1*A

% The third Householder matrix to the third column

v=A2(3:4,3);
v(1)=v(1)+norm(v);
H3=blkdiag(1,1,eye(2)-2*v*v'/(v'*v)) 

% this zeros out the third column below the main diagonal

R=H3*H2*H1*A % This is already an upper triangular matrix that can be denoted by R.

% The orthogonal matrix Q can be constructed as follows.

Q=H1*H2*H3 % Because this matrices are symmetric.

% Compare these Q and R matrices with the result of Example 1. We can see
% that these matrices are very similar to those (some signs are different only)
% and similarly A=QR (the QR decomposition is not unique in general).

%% Givens rotations
% There is another possibility to transform a matrix to upper triangular form with orthogonal matrices. We can use rotations instead of reflections. The matrix (orthogonal !) of the rotation with angle \theta can be seen below. If the angle \theta is chosen such that the second element of the image vector is zero then we say that we have defined a Givens rotation.
%
%
% Example 5. (Givens rotation)

% Define the Givens rotation matrix to the vector x=[2;5] and to the last
% two elements of the vector y=[1;3;-5;2;5]

x=[2;5]
y=[1;3;-5;2;5]

% First we compute the sin and cos values of the rotation angle (we do not need the 
% angle explicitly)

c=x(1)/norm(x), s=-x(2)/norm(x)

% and construct the Givens rotation matrix

G=[c,-s;s,c]

% The effect of G is

G*x 

% multiplication with G zeros out the second element of x (the 2-norm does not 
% change)

% For the vector y we can use the same 2D rotation as follows

blkdiag(1,1,1,G)*y

% this transformation zeros out the last element of the vector y.

% Because the Givens rotation matrices are orthogonal and they introduce zeros into a vector, they can be used to construct the QR decomposition of a matrix (this is the second method of the production of the QR decomposition). We have to apply the technique of Example 5 in the following way.
%
% Example 6. (QR decomposition with Givens rotations)

% Let us give the QR decomposition of the matrix A=[1,4;2,6;0,-2;3,7] with 
% Givens rotations

A=[1,4;2,6;0,-2;3,7]

% Givens rotation to the last two elements of the first column.

x=A(3:4,1);
c=x(1)/norm(x); s=-x(2)/norm(x);

G1=blkdiag(1,1,[c,-s;s,c])

% The result of the first transformation

A1=G1*A % the last element in the first column is zero

% Givens rotation to the 2. and 3. elements of the first column.

x=A1(2:3,1);
c=x(1)/norm(x); s=-x(2)/norm(x);

G2=blkdiag(1,[c,-s;s,c],1)

% The result of the second transformation

A2=G2*G1*A % the last two elements in the first column are zero

% Givens rotation to the 1. and 2. elements of the first column.

x=A2(1:2,1);
c=x(1)/norm(x); s=-x(2)/norm(x);

G3=blkdiag([c,-s;s,c],1,1)

% The result of the third transformation

A3=G3*G2*G1*A % the last three elements in the first column are zero

% Now comes the column to zero out

% Givens rotation to the last two elements of the second column.

x=A3(3:4,2);
c=x(1)/norm(x); s=-x(2)/norm(x);

G4=blkdiag(1,1,[c,-s;s,c])

% The result of the fourth transformation

A4=G4*G3*G2*G1*A % the last element in the second column is zero, the first 
% column does not change

% Givens rotation to the 2. and 3. elements of the second column.

x=A4(2:3,2);
c=x(1)/norm(x); s=-x(2)/norm(x);

G5=blkdiag(1,[c,-s;s,c],1)

% The result of the fifth transformation

R=G5*G4*G3*G2*G1*A % this is already an upper triangular matrix, this is why we 
% denoted it by R.

% The Q matrix can be obtained as follows.

Q=G1'*G2'*G3'*G4'*G5'

% This is a QR decomposition indeed

Q*R % gives back A

%
% We can see from the number of operations that the CPU time of a QR decomposition for a full matrix with Givens rotations is 50% larger than that of the Householder QR. This is why we generally use Householder's method to produce the QR decomposition. Givens' method is efficient e.g. for upper Hessenberg matrices.
%
%
%
%% Over-determined linear systems
%
% Now we consider linear systems with matrices with full column rank. Generally, these are linear systems where we have more equations than unknowns.
%
% Thus the x_{\mathrm{LS}} solution (the subscript here means solution in the "least square" sense) minimizes the 2-norm of the residual r=b-\mathrm{Ax}. The x_{\mathrm{LS}} solution can be obtained with two methods: solving the so-called normal equation or using the QR decomposition of the coefficient matrix.
% Solution with the normal equation
%
% Solution with the QR decomposition
% Here we use that multiplications with orthogonal matrices do not change to 2-norm.
%
%
% Example 7. (Solution of over-determined linear systems)

% Solve the over-determined system with A=[1,4;2,6;0,-2;3,7] and b=[1;2;3;4].

A=[1,4;2,6;0,-2;3,7] 
b=[1;2;3;4]

% Normal equation

x=(A'*A)\(A'*b) % This is the x_LS solution


% QR decompotition

[Q,R]=qr(A) % First we compute the QR decomposition of the matrix A. Here we could use 
% the result of Example 6.

d=Q'*b;
d=d(1:2)
x=R(1:2,1:2)\d % We have to take the upper square part of the matrix R (2x2 block) and 
% the first 2 element of the vector Q^Tb. This gives x_LS again.

% Solution with Matlab's built-in function
% The left division function can solve over-determined systems too.

x=A\b % This gives x_LS again.

%
%% Further problems
%
% Solve problems 38-42. from our problem booklet. Try to solve the problems both in Matlab (using e.g. the previous codes) and manually. In problem 39 use also Householder reflections to obtain the QR decomposition. The solutions can be found below (and also in the m-file of the course) but first give it a try.
%
% Problem 38

x=[2;6;-3]; % we define the vector

% Householder reflection
v=x; v(1)=v(1)+norm(v) % the normal vector of the reflection plane
H=eye(3)-2*v*v'/(v'*v) % the Householder matrix
H*x % Matrix H zeros out the elements with index greater than 1.

% Givens rotation
c=6/norm(x(2:3)); s=3/norm(x(2:3));
G=blkdiag(1,[c,-s;s,c]) % The Givens rotation matrix
G*x % Matrix G zeros out the last element of the vector.


% The manual computations should give the following results (obtained by the symbolic toolbox).

x=sym([2;6;-3]); % we define the vector

% Householder reflection
v=x; v(1)=v(1)+norm(v) % the normal vector of the reflection plane
H=eye(3)-2*v*v'/(v'*v) % the Householder matrix
H*x % Matrix H zeros out the elements with index greater than 1.

% Givens rotation
c=6/norm(x(2:3)); s=3/norm(x(2:3));
G=blkdiag(1,[c,-s;s,c]) % The Givens rotation matrix
G*x % Matrix G zeros out the last element of the vector.

%
% Problem 39.

A=[0,0;1,3;0,2]

% Givens rotation

x=A(1:2,1);
c=x(1)/norm(x), s=-x(2)/norm(x)
G1=blkdiag([c,-s;s,c],1) % The first Givens matrix
A1=G1*A
x=A1(2:3,2);
c=x(1)/norm(x), s=-x(2)/norm(x)
G2=blkdiag(1,[c,-s;s,c]) % The second Givens matrix
R=G2*A1 % R
Q=G1'*G2' % Q
Q*R % Check whether A=QR 

% Householder reflection


v=A(:,1);
v(1)=v(1)+norm(v)
H1=eye(3)-2*v*v'/(v'*v)

% this zeros out the first column of A below the diagonal

A1=H1*A

% The second Householder matrix to the second column

v=A1(2:3,2);
v(1)=v(1)+norm(v);
H2=blkdiag(1,eye(2)-2*v*v'/(v'*v)) 
% this zeros out the second column below the main diagonal

R=H2*H1*A % This is already R.

Q=H1*H2 % This is matrix R.

% Built-in function

[Q,R]=qr(A) % Computation of Q and R with the built-in 
% function, these matrices are the same as we obtained in the previous two parts.

% Here there is no difference between Matlab and manual solutions.

%
% Problem 40
% Assume that there are two different A=\mathrm{QR}=Q_1 R_1 decompositions with positive diagonals in R and R_1. Then we have
% Q^{-1} Q_1 =RR_1^{-1} (let us denote this matrix by D).
% Here the left hand side is an orthogonal matrix. Its inverse is its transpose. The right hand side is an upper triangular matrix. Its inverse is an upper triangular matrix again. This is possible only if D is a diagonal matrix. At the same time, this matrix must be orhogonal, which means that all its diagonal elements must be +1 or -1.
% Thus, R={D\;R}_1. Because all diagonal elements in R and in R_1 are positive, D is the identity matrix, which gives the contradiction R=R_1.
%
% Problem 41

A=[0,0;1,3;0,2]  % the matrix from problem 39.
b=[3;4;1];

% built in solver

A\b

% normal equation

(A'*A)\(A'*b)

% QR

[Q,R]=qr(A);
c=Q'*b;
R(1:2,1:2)\c(1:2)

% Here there is no difference between Matlab and manual solutions.

%
% Problem 42
% The value to be minimized can be written in the form
% {\left\|\left\lbrack \begin{array}{ccc}
% 1 & 1 & 1\\
% 0 & 0 & 1\\
% 1 & 1 & 1\\
% 4 & 2 & 1
% \end{array}\right\rbrack \left\lbrack \begin{array}{c}
% a\\
% b\\
% c
% \end{array}\right\rbrack -\left\lbrack \begin{array}{c}
% 1\\
% 0\\
% 3\\
% 2
% \end{array}\right\rbrack \right\|}_2^2, which exactly means that we have to compute the x_LS solution of the linear system with matrix \left\lbrack \begin{array}{ccc}
% 1 & 1 & 1\\
% 0 & 0 & 1\\
% 1 & 1 & 1\\
% 4 & 2 & 1
% \end{array}\right\rbrack and right hand side \left\lbrack \begin{array}{c}
% 1\\
% 0\\
% 3\\
% 2
% \end{array}\right\rbrack.

A=[1,1,1;
   0,0,1;
   1,1,1;
   4,2,1]

b=[1;0;3;2]

x=A\b  % We can use here the built-in solver.

% Thus the best polynomial is p(x)=-x^2+3x

% Just to see the points and the fitted polynomial

figure
fplot(@(x)-x.^2+3*x,[0,2])
hold on
plot([1,0,1,2],[1,0,3,1],'ro') 
hold off
snapnow

%
%% Homework - remember to solve the homework from the webpage! The deadline is 22 March.
%
%