%% Esercizi sulle scomposizioni matriciali

% Generare una matrice di numeri casuali denominata X di dimensione nxp (ad esempio
% n=100 e p=5). 
% Calcolare il rango della matrice X tramite la funzione rank
% v. file 

% Calcolare la matrice di covarianze S.
% Verificare la scomposizione spettrale della matrice 
%  S=V \Lambda V'.
% Verificare l'ortogonalità della matrice V
% Verificare che la somma degli autovalori coincide con la traccia della
% matrice S.
% Verificare che il prodotto degli autovalori coincide con il determinante della
% matrice S.

% Riordinare le colonne della matrice V e La (\Lambda) in modo tale che
% La(1,1) sia il grande autovalore e V(:,1) sia l'autovettore associato
% La(2,2) sia il secondo più grande autovalore e V(:,2) sia l'autovettore associato
% ......

% Verificare tramite la funzione rank che la matrice pxp 
% v_i v_i^' abbia rango 1, i=1, 2, ..., p 
% Ricostruire la matrice S tramite la somma di p matrice quadrate ciascuna
% di rango 1
% $$\sum_i=1^p \lambda_i v_i v_i^'$$

% Effettuare la scomposizione in valori singolari della matrice Xtilde
% Xtilde= U \Gamma Vtramitesvd^'
% Verificare che la matrice Vtramitesvd è uguale (a meno del segno) alla matrice V ottenuta tramite la
% scomposizione spettrale della matrice S
% Verificare che (\Gamma/sqrt(n-1)).^2 = \Lambda dove \Lambda è la matrice
% degli autovalori ottenuta in precendenza

% OSSERVAZIONE: tenere presente che nella scomposizione spettrale gli autovalori sono ordinati
% dal più piccolo al più grande
% nella scomposizione in valori singolari i valori singolari (radici
% quadrate degli autovalori) sono ordinati dal più grande al più piccolo

% Controllare che la matrice di dimensione nxp $\gamma_i u_i v_i^'$ 
% abbia rango 1, i=1, 2, ..., p

% Ricostruire la matrice Xtilde tramite la somma di p matrici di dimensioni nxp
% ciascuna di rango 1 
% $\sum_i=1^p \gamma_i u_i v_i^'$

% Calcolare la somma dei quadrati delle differenze che si ottengono quando
% la matrice Xtilde viene approssimata dalla matrice Xcappello
% Xcappello=$\sum_i=1^2 \gamma_i u_i v_i^'$
% ossia dalla matrice di rango 2 contenente i due autovalori più grandi 
% Verificare che la somma dei quadrati delle differenze è uguale alla somma
% dei quadrati dei valori singolari (somma che va da 3 a p).
% Calcolare la somma dei quadrati delle differenze in maniera matriciale e
% non matriciale 
% Osservazione: data una matrice X la somma dei quadrati degli elementi
% di X in maniera matriciale si scrive = traccia(X' X)=traccia (X X')


%% Scomposizione in valori singolari di una matrice quadrata e simmetrica
n=100;
p=5;
X=randn(100,5);

disp('Rango della matrice X')
disp(rank(X))

% Calcolare la matrice di covarianze S.
S=cov(X);
[V,La]=eig(S);

% Verificare la scomposizione spettrale della matrice 
%  S=V \Lambda V'.
% S ricostruita tramita la scomposizione spettrale
Schk=V*La*V';


%% Controllo dell'ortogonalità della matrice V
disp(V'*V)

%% La somma degli autovalori coincide con la traccia
disp('Verifico che la traccia della matrice S è uguale alla somma dei suoi autovalori')
disp(['Traccia della matrice S=' num2str(trace(S))])
disp(['Somma degli autovalori della matrice S=' num2str(sum(diag(La)))])

%% Il prodotto degli autovalori coincide con il determinante
detS=det(S);
proLa=prod(diag(La));
disp('Verifico che il prodotto degli autovalori della matrice S coincide con il determinante')
disp(['Traccia della matrice S=' num2str(detS)])
disp(['Prodotto degli autovalori della matrice S=' num2str(proLa)])


%% Riordinare la matrici V e Lambda 
la=diag(La);
[~,indsor]=sort(la,'descend');
% Riordino le colonne della matrice V e La in modo tale che
% La(1,1) sia il grande autovalore e V(:,1) sia l'autovettore associato
% La(2,1) sia il secondo più grande autovalore e V(:,2) sia l'autovettore associato
V=V(:,indsor);
lasor=la(indsor);
La=diag(lasor);

% Controllo il rango della matrice
% v_i v_i^' quando i=2
v2=V(:,2);
disp('Rango della matrice v_2 v_2^')
disp(rank(v2*v2'))


Schk1=zeros(p,p);
% S ricostruita tramite la somma di p matrici
% quadrate ciascuna di rango 1
for j=1:p
    Schk1=Schk1+La(j,j)*V(:,j)*V(:,j)';
end

%% Scomposizione in valori singolari della matrice Xtilde
% Xtilde = matrice degli scostamenti dalla media)
Xtilde=X-mean(X);

[U,Gam,Vtramitesvd]=svd(Xtilde,'econ');
Xtildechk=U*Gam*Vtramitesvd';
disp('Ricostruzione della matrice Xtilde')
disp('Controllo che il massimo delle differenze in valore assoluto sia un numero piccolissimo')
max(max(abs(Xtilde-Xtildechk)))


Xtildechk1=zeros(n,p);

for j=1:3
    Xtildechk1=Xtildechk1+U(:,j)*Gam(j,j)*Vtramitesvd(:,j)';
end
disp('Controllo che il massimo delle differenze in valore assoluto sia un numero piccolissimo')
max(max(abs(Xtilde-Xtildechk1)))



% Verifica che (\Gamma/sqrt(n-1)).^2 = \Lambda
disp('(\Gamma/sqrt{ n-1})^2')
disp((Gam/sqrt(n-1)).^2)
disp('Matrice \Lambda degli autovalori della matrice di covarianze')
disp(La)

% OSSERVAZIONE: nella scomposizione spettrale gli autovalori sono ordinati
% dal più piccolo al più grande
% nella scomposizione in valori singolari i valori singolari (radici
% quadrate degli autovalori) sono ordinati dal più grande al più piccolo
% ecco perché  è stato necessario precedentemente effettuare il
% riordinamento degli autovalori

disp('Rango della matrice di dimensione nxp u_3 v_3^''')
u3=U(:,3);
v3=V(:,3); % oppure v3=Vtramitesvd(:,3)
disp(rank(u3*v3'))



%% Calcolare la somma dei quadrati delle differenze che si ottengono quando
% la matrice Xtilde viene approssimata dalla matrice Xcappello
% Xcappello=$\sum_i=1^2 \gamma_i u_i v_i^'$
% ossia dalla matrice di rango 2 contenente i due autovalori più grandi 

% Calcolo  Xcappello=\sum_i=1^2 \gamma_i u_i v_i^'
Xcappello=zeros(n,p);
for j=1:2
    Xcappello=Xcappello+U(:,j)*Gam(j,j)*Vtramitesvd(:,j)';
end
% Calcolo la somma dei quadrati delle differenze tra la matrice Xtilde e la
% matrice Xcappello (in maniera non matriciale)
disp(sum(sum((Xtilde-Xcappello).^2)))

% Calcolo la somma dei quadrati delle differenze tra la matrice Xtilde e la
% matrice Xcappello (in maniera matriciale tramite la traccia)
disp(trace((Xtilde-Xcappello)'*(Xtilde-Xcappello)))

% Verifico che la somma dei quadrati delle differenze è uguale alla  somma
% dei quadrati dei valori singolari da 3 a p
disp(sum(diag(Gam(3:end,3:end)).^2))

% La matrice Xcappello è la miglior rappresentazione di rango 2 della
% matrice Xtilde di partenza nel senso che se si sostituisce alla matrice
% Xcappello qualsiasi altra matrice di rango 2, la somma dei quadrati
% risulta maggiore.