%% Time Varying Systems and Computation
% Klaus Diepold, Sven Gronauer, Simon Schumann
%
%% 1. Linear Algebra
%% 1.1 Warm-up
% How can we construct data sparse and structured matrices?
% diagonal matrix:

d=diag([1,2,3,4])

% triangular matrix:

diag([1,2,3])+diag([4,5],-1)+diag(6,-2)

% sparse matrix:

sprand(5,5,0.2) %random 500-by-1000 sparse matrix with density 0.2

% Hankel matrix:

hankel([1,2,3,4],[4,5,6,7])

% Toeplitz matrix:

toeplitz([1,4,3,2],[1,2,3,4])

% Vandermonde matrix:

fliplr(vander([1 2 3 4]))

%
% Let's look again at our diagonal matrix d. Currently, all 16 elements of the matrix are stored, which is not very efficient considering that most of the entries are zero. To fix that, we can use the sparse dataformat provided by Matlab, such that zeros are not stored explicitly anymore. Instead, only the nonzero elements are stored, together with the information about their position in the matrix. So let's save some memory while using sparse matrices:

sparse(d)

%
%% 1.2 Cost of linear algebra
% Let's plot the time consumed for computing matrix-vector product for dense square matrices of different sizes.

dimensions = 1000:1000:20000;
dense_times = zeros(size(dimensions));
time_index = 1;
for n=dimensions
    A = randn(n, n);
    b = randn(n, 1);
    dense_times(time_index) = timeit(@()A*b);
    time_index = time_index + 1;
end
plot(dimensions, dense_times);

% What can you tell about the asymptotic complexity of computing matrix-vector product for dense square matrices by looking at the plot?
% Hint: Add linear and quadratic trend lines to the plot by using "plotLinearTrend(dimensions,dense_times)" and "plotQuadraticTrend(dimensions,dense_times)".
% Let's repeate this for sparse square matrices of different sizes.  (Each matrix has 100.000 nonzero elements).

sparse_times = zeros(size(dimensions));
time_index = 1;
for n=dimensions
    density = 100000 / (n^2);
    A = sprand(n, n, density);
    b = randn(n, 1);
    sparse_times(time_index) = timeit(@()A*b);
    time_index = time_index + 1;
end
plot(dimensions, sparse_times);

% What can you tell about its asymptotic complexity?
%
%% 1.3 Solving systems of linear equations
% Consider the following system of equations

A=randn(5,3)*randn(3,4);
x = randn(4,1);
b= A*x;

% What can we tell about the matrix \mathbf{A} by looking at how we construct it?
% To solve this equation with respect to x, we multiply \mathbf{b} with \mathbf{A}^{\dagger} = (\mathbf{A}^{T}\mathbf{A})^{-1}\mathbf{A}^{T}, where \mathbf{A}^{\dagger}  is the pseudo inverse of \mathbf{A}.

A_pinv = inv(A'*A)*A';
x_p = A_pinv * b;

% To evaluate how well \mathbf{x} is predicted by \mathbf{x}_{p}, let's inspect the mean squared error of the prediction.

fprintf("MSE: %f", mean((x_p - x).^2))

% Run the code several times. What do you observe?
% Instead of calculating pseudo inverse by ourselves let's use the matlab command pinv to solve this linear equation.

A_pinv = pinv(A);
x_p = A_pinv * b;

%  Let's again inspect the mean squared error of the prediction.

fprintf("MSE: %f", mean((x_p - x).^2))

% Do you see any differences? If yes, what could be the reason for this difference?
%
%% 1.4 Questions
% How can you save memory while using sparse matrices?
% What can you tell about the asymptotic complexity of computing matrix-vector product for dense square matrices by looking at the plot?
% What can you tell about the asymptotic complexity of computing matrix-vector product for sparse square matrices by looking at the plot?
% What can you tell about the matrix \mathbf{A} at line 29 by looking at how we construct it?
% Run the code several times. What do you you observe about the mean squared error computed at line 34?
%  What is the reason for the difference between the two mean squared errors computed at lines 34 and 37?