%% Time-Varying Systems and Computations
% Klaus Diepold, Matthias Kissel, Sven Gronauer, Nehir Güzelkaya
%
%% 1. Separable Matrix
% We can create a separable matrix by computing outer product of two random vectors p and q:

n = 6;
p = randn(n, 1);
q = randn(n, 1);
separable_mat=p*q'

% Now, let's check the rank of the separable matrix:

rank(separable_mat)

% What do you observe? Is this the case for all separable matrices?
%
%% 2. Semiseparable Matrix
% We can create a random semiseparable matrix by combining two rank 1 matrices:

n = 6;
p = randn(n, 1);
q = randn(n, 1);
v = randn(n, 1);
w = randn(n, 1);
semiseparable_mat = tril(p*q') + triu(v*w', 1)

% Let's again check its rank:

rank(semiseparable_mat)

% Is it different than the rank of a separable matrix?
% Finally, let's check if all submatrices taken out of the strictly upper or lower diagonal matrix are of low rank (i.e. their rank is <= 1 in our case):

lower_ranks = zeros(1, n-1);
upper_ranks = zeros(1, n-1);
for i=1:n-1
    lower_block = semiseparable_mat((i+1):n, 1:i);
    lower_ranks(i) = rank(lower_block, 1e-3);
    upper_block = semiseparable_mat(1:i, (i+1):n);
    upper_ranks(i) = rank(upper_block, 1e-3);
end
result = all(lower_ranks <= 1) && all(upper_ranks <= 1)

% At the end of the script, this part of the code is extended into a function are_submatrices_low_rank() that we are going to use in the following sections.
%
%% 3. Matrix-Vector Multiplication with Semi-Separable Matrix
% We now have a look at a matrix-vector multiplication, where the coefficient matrix is a lower triangular semi-separable matrix.

n = 5;
p = randn(n, 1);
q = randn(n, 1);
lt_semiseparable_mat = tril(p*q');

% Now, we can compute the matrix vector mutliplication y of lt_semiseparable_mat with a random vector u as

u=randn(n,1);
y=zeros(n,1);
sum=0;
for i=1:n
    sum=sum+q(i)*u(i);
    y(i)=p(i)*sum;
end

% Finally, we can check if our result is true:

y
y_ref=lt_semiseparable_mat*u
mse=norm(y-y_ref)^2

%
%%  4. Inverse of Bi-Diagonal and Tri-Diagonal Matrices
% We first create Bi-Diagonal and Tri-Diagonal Matrices:

n=5;
bidiag_mat = diag(randn(n, 1)) + diag(randn(n-1, 1), -1);
display(bidiag_mat);
tridiag_mat = diag(randn(n, 1)) + diag(randn(n-1, 1), -1) + diag(randn(n-1, 1), 1);
display(tridiag_mat);

% Now, let's check if the submatrices taken out of Bi-Diagonal and Tri-Diagonal Matrices are of low rank.

display(are_submatrices_low_rank(bidiag_mat));
display(are_submatrices_low_rank(tridiag_mat));

% Finally, we can check if this also applies to their inverses.

inverse_bidiag_mat = bidiag_mat^-1;
display(are_submatrices_low_rank(inverse_bidiag_mat));
inverse_tridiag_mat = tridiag_mat^-1;
display(are_submatrices_low_rank(inverse_tridiag_mat));

%
%% 5. QR Decomposition of a semiseparable matrix
% Let's compute the QR Decomposition of the semiseparable matrix from section 1.

[ssQ, ssR] = qr(semiseparable_mat);

% Now, let's check if the submatrices taken out of the resulting Q matrix are still of low rank.

display(are_submatrices_low_rank(ssQ));

% Are the submatrices taken out of the resulting R matrix also of low rank?
%
% Function that checks if all submatrices taken out of the strictly upper or lower diagonal matrix are of low rank i.e. their rank is <= 1:

function result = are_submatrices_low_rank(mat)
    n = size(mat,1);
    lower_ranks = zeros(1, n-1);
    upper_ranks = zeros(1, n-1);
    for i=1:n-1
        lower_block = mat((i+1):n, 1:i);
        lower_ranks(i) = rank(lower_block, 1e-3);
        upper_block = mat(1:i, (i+1):n);
        upper_ranks(i) = rank(upper_block, 1e-3);
    end
    result = all(lower_ranks <= 1) && all(upper_ranks <= 1); 
end