%% Tutorial 2. Exploring Matrix Inversion when Solving Linear Systems 
%
% In this tutorial, we will see why performing matrix inversion is generally 
% not a good idea. A few simulation experiments will be conducted to
% demonstrate the accuracy and speed of two solution methods for a linear
% system of equations. 
%

%% 1. Direct Solution vs. Matrix Inversion
% In general, matrix inversion arises when trying to solve the following
% linear system A*x = b via x = A^-1 * b
% However, the computation of the inverse can result in long execution 
% time and poor numerical accuracy, so the matrix division operator (or
% backslash) is preferred. The backslash operator performs Gaussian
% elimination of the system, without forming the inverse.  
%
% In this section, we will carry out the solution of a linear system by way
% of these two methods and compare the results of the calculation. 

% Specify vector of matrix dimensions to use for the calculation. 
ndim = round(2.^(7:.3:11)); % Note that this gives you exponential spacing.
reps = 10; % Repeat the calculation 10 times. 

% We start by creating arrays that will store computation time and error of
% the calculation. 
t = zeros(length(ndim),2);
err = zeros(length(ndim),2);

% Write a loop that performs the computation for linear systems of
% different dimensions, and computes the error of 10 different trials. 
for k = 1:length(ndim)
    
    err(k,:) = 0;
    t(k,:) = 0;
    
    for r = 1 : reps
        A = randn(ndim(k));
        y = randn(ndim(k),1);
        b = A*y;
        tic
        x = inv(A)*b;
        t(k,1) = t(k,1) + toc;
        err(k,1) = err(k,1) + norm(x-y)/ndim(k)/norm(y);
        tic
        x = A\b;
        t(k,2) = t(k,2) + toc;
        err(k,2) = err(k,2) + norm(x-y)/ndim(k)/norm(y);
    end
    
    err(k,:) = err(k,:)/reps;
    t(k,:) = t(k,:)/reps;
    
end


% Display the results of the calculation in a subplot. 
figure(1);clf;
subplot(211);
loglog(ndim,t);
xlabel('Matrix Size n');
ylabel('Time(sec)');
legend('Matrix Inversion','Direct Solution',-1)
grid on
ax = axis;
ax(2) = 4000;
axis(ax);

subplot(212);
loglog(ndim,err);
xlabel('Matrix Size n');
ylabel('||x-y||');
legend('Matrix Inversion','Direct Solution',-1)
grid on;
ax = axis;
ax(2) = 4000;
axis(ax);

%% 2. Data Fitting (Overdetermined Equations)
% Another situation, which we will study in greater detail later in the
% course is fitting of data. Consider a problem where you have more
% experimental data points than prediction variables.  In this case you
% can't solve it exactly, but instead solve it in a "least square" sense.
% Note that the fitted coefficients don't exactly match up with the values
% we use to generate the data.

% Suppose we want to fit a line to noisy data:
t = [0:.1:1]';
A = [t,ones(size(t))];
b = 3*t -1 +0.2*randn(size(t));

% Now, solve the system. 
x = A\b

figure(2)
plot(t,b,'o',t,A*x)
xlabel('t')
ylabel('b')
title('Fitting a line to data')

%% 3. Underdetermined Equations
% Another case will be when we have more variables to predict than data to
% support it. In this case there are an infinite number of solutions.  We
% get different answers depending on the algorithm you use.
A = magic(8);
A = A(1:6,:);
b = sum(A,2);

% Check the dimensions of A.  
sizeA = size(A)
rankA = rank(A)

% We will get a warning in this case because of the rank of the matrix.
% Alternatively, we can use the command pinv in MATLAB, which generates a
% square matrix and then solves the system. The length of the solution
% vector is minimized using pinv.

xbackslash = A\b
 

xpinv = pinv(A)*b
 

%% 4. Simple Example of an Ill-Conditioned Matrix
% Recall from your linear algebra course that some matrices may be
% ill-conditioned if a small change in their constant coefficients results
% in a large change in their solution. Let us analyze one such matrix.
% Because the condition number is very large, the solutions will be
% significantly different from the case with no error!

A = [11  10  14;
    12  11 -13;
    14  13 -66];

b = [1;1;1];

% We have the condition number:
condA = cond(A)

% We introduce errors in A and b:
Ap = A;
Ap(2,2) = Ap(2,2) + 0.0001;
bp = [1.0001; 0.9999; 1.0001];

% We now see how the answer changes: 
x= A\b
 
% Error in b:
xberror = A\bp

% Error in A
xAerror = Ap\b