%%
%% Time Varying Systems and Computation
% Klaus Diepold, Sven Gronauer, Simon Schumann
%
%% 4. Matrix decompositions
%% 4.1 Tools
% Givens Rotation:  The Givens elementary transformation is:
%
% It represents an elementary rotation over an angle \phi:
%
% We can choose the value for this parameter \phi such that the resulting transformation annihilates a predetermined element in a row (with c\;\dot{=} \;\cos \;\phi and s\;\dot{=} \;\sin \;\phi)
% \left\lbrack \begin{array}{cc}
% c & -s\\
% s & c
% \end{array}\right\rbrack \left\lbrack \begin{array}{cccc}
% a_1  & a_2  & \ldotp \ldotp \ldotp  & a_m \\
% b_1  & b_2  & \ldotp \ldotp \ldotp  & b_m
% \end{array}\right\rbrack =\left\lbrack \begin{array}{ccc}
% {c\cdot \;a}_1 -{s\cdot b}_1  & {c\cdot \;a}_2 -{s\cdot b}_2  & \ldotp \ldotp \ldotp \\
% s\cdot a_1 +c\cdot b_1  & s\cdot a_2 +c\cdot b_2  & \ldotp \ldotp \ldotp
% \end{array}\right\rbrack \dot{=} \;\left\lbrack \begin{array}{cccc}
% \sqrt{\left|a_1 {\left|\right.}^2 +\right|b_1 {\left|\right.}^2 } & \star  & \ldotp \ldotp \ldotp  & \star \\
% 0 & \star  & \ldotp \ldotp \ldotp  & \star
% \end{array}\right\rbrack
% when s\cdot a_1 +c\cdot b_1 =0 or \tan \;\phi =-\frac{b_1 }{a_1 }.
% Code example:

A = randn(2,2)
phi = -atan(A(2,1) / A(1,1));
R = [cos(phi), -sin(phi); sin(phi), cos(phi)];
R * A

% Check the determinant of the Givens rotation matrix

det(R)

% Is this the case for all rotation matrices?
% Householder Reflection: Householder reflections can be used to reflect a vector onto a standard basis vector.
%
% In the picture, the vector a is reflected onto the basis vector e_1. The resulting vector b has the same length as a.
% Code example in the three-dimensional space:

a = randn(3,1);
e1 = [1;0;0];

% The vector b is obtained by scaling e_1 ={\left\lbrack 1\;0\;0\right\rbrack }^T with the norm of a.

b = e1 * norm(a);

% Now we find v:=2w as  the difference between a and b. Note that we first make sure that the sign of b_1 is different from the sign of a_1 (imagine what otherwise happens if a=b).

b(1) = b(1) * ((a(1) < 0)*2 - 1);
v = a - b;

% Then, we calculate the vector w by projecting a onto the axis spanned by u. Hereby u is a vector of length 1 with the same direction as v.

u = v / norm(v);
P = u*u';
w = P*a;

% Finally, we subtract 2w from a to obtain b.

b_proj = a - 2*w
b - b_proj
norm(a)
norm(b_proj)

% A sequence of Householder transformations can be used to introduce subdiagonal zeros into a matrix ---> Homework
%
%% 4.2 QR decomposition
% Consider the following matricies

A_square = randn(4,4);
A_tall = randn(5,3);
A_wide = randn(3,5);

% and their QR decompositions, computed with the Matlab command qr():

[Q_square, R_square] = qr(A_square)
[Q_tall, R_tall] = qr(A_tall)
[Q_wide, R_wide] = qr(A_wide) 

% Lets check out the structure of the 'R' matrices and the 'Q' matrices properties.

figure
subplot(2,2,1)
spy(R_square)
title("R square")
subplot(2,2,2)
spy(R_tall)
title("R tall")
subplot(2,2,[3, 4])
spy(R_wide)
title("R wide")

% Square matrix Q
disp("Results for square matrix")
Q_square * Q_square'
Q_square'* Q_square
% Tall matrix Q
disp("Results for tall matrix")
Q_tall * Q_tall'
Q_tall'* Q_tall
% Wide matrix Q
disp("Results for wide matrix")
Q_wide * Q_wide'
Q_wide'* Q_wide

% Now let's try what happens if we set the second parameter of the qr() to 'econ'

[Q_square, R_square] = qr(A_square,'econ')
[Q_tall, R_tall] = qr(A_tall,'econ')
[Q_wide, R_wide] = qr(A_wide,'econ') 

% And finally have a look into the structure of the new 'R' matrices and  'Q' matrices properties.

figure
subplot(2,2,1)
spy(R_square)
title("R square")
subplot(2,2,2)
spy(R_tall)
title("R tall")
subplot(2,2,[3, 4])
spy(R_wide)
title("R wide")

% Square matrix Q
disp("Results for square matrix")
Q_square * Q_square'
Q_square'* Q_square
% Tall matrix Q
disp("Results for tall matrix")
Q_tall * Q_tall'
Q_tall'* Q_tall
% Wide matrix Q
disp("Results for wide matrix")
Q_wide * Q_wide'
Q_wide'* Q_wide

% Do you see any differences?
%% 4.3 QR decomposition using Givens rotation
% Embedded Givens Rotations can be utilized to compute the QR decomposition of a matrix. The procedure starts at the first column and last row. It anhilates all elements except the first one in the first column and then continues with the next column where all elements except the first and second one get anhilated (starting from the bottom row), and so on.
%
% The end result is:
%
% Code example:

A = randn(5,5)     % Matrix which we want to factorize
Q = eye(size(A,1)) % Initial value for Q factor
R = A              % Initial value for R factor
figure
for i=1:size(A,2)  % Loop over columns
    for j=size(A,1):-1:i+1 % Loop over rows (start from bottom)
        phi = -atan(R(j,i) / R(j-1,i));
        G_embedded = eye(size(R,1)); % Embedded givens rotation for j'th row
        G_embedded([j-1,j],[j-1,j]) = [cos(phi), -sin(phi); sin(phi), cos(phi)];
        R = G_embedded * R;  % Anhilate the next element, producing R
        Q = Q * G_embedded'; % "Record" transformations used in fator Q
        subplot(size(A,1),size(A,2), (i-1)*size(A,2) + (size(A,1) - j + 1))
        spy(abs(R) > 1e-10)
    end
end

% QR and QL decompositions using Householder transformations, LQ and RQ decompositions using Givens rotations ---> Homework
%% 4.4 Questions
% What can you say about the determinant of rotation matrices?
% Is there any consequences of using economic qr decomposition in the orthogonality of Q matrices?