% lecture 16, question 2.

clc;
clear;
close all;
warning off;

for loop=1:5
  % set color map
  colour = ['r';'g';'b';'k';'m'] ;

  M = 5;

  A = randn(M);
  [U,S,V] = svd(A); 
  clear A S;

  s = rand(M,1);
  s = sort(s);
  s = flipud(s);
% (c)
%  s = s.^6;
  S = diag(s);
  clear s;

  A = U*S*V';

  [U2,S2,V2] = svd(A);
  scale = sign(diag(U2'*U));
  U2 = U2*diag(scale);
  V2 = V2*diag(scale);
  
  normU = norm(U - U2);
  normS = norm(S - S2);
  normV = norm(V - V2);
  normA = norm(A - U2*S2*V2');

  disp(['Example #' num2str(loop)]);  
  disp(['The 2-Norm of U is ' num2str(normU)]); 
  disp(['The 2-Norm of S is ' num2str(normS)]);  
  disp(['The 2-Norm of V is ' num2str(normV)]);  
  disp(['The 2-Norm of A is ' num2str(normA)]);  
  disp('*********************************');

  figure(1);
  plot(1:M',diag(U2'*U),colour(loop));
  hold on
  title('T = diag(U2''*U)');
  xlabel('j')
  ylabel('T(j)');

  figure(2);
  plot(1:M',diag(V2'*V),colour(loop));
  hold on
  title('R = diag(V2''*V)');
  xlabel('j')
  ylabel('R(j)');

end