% Simple code to demonstrate Theorem 9.1 in 
% Halko/Martinsson/Tropp (page 267)

 m = 200;
 n = 50;

 kmax = 20;
 dat = zeros(kmax,2);

 U = orth(randn(m,n));             % random left singular vectors
 V = orth(randn(n));               % random right singular vectors
 Sig = diag(2.^[0:-1:1-n]);        % exponentially decaying singular values
 A = U*Sig*V';                     % construct A with specified singular values
 
 for k=1:kmax                      % iterate over reduction dimension, k
    ell = 2*k;                     % enlarged space

    Omega = randn(n,ell);          % random sample matrix

    Y = A*Omega;                   % action of A on random samples

    At = U'*A;                   
    Yt = U'*Y;                     % equivalent projected coordinates

    Q = orth(Yt);
    PiYt = Q*Q';                   % orthogonal projector onto R(Yt)

    Im = eye(m);
    In = eye(n);

%    fprintf('True Error: norm((In-PiYt)*A) = %20.15f\n', norm((In-PiYt)*At))
    dat(k,1)= norm((In-PiYt)*At);  % exact error, = norm((I-PiY)*A)

    Sig2 = Sig(k+1:end,k+1:end);   % neglected singular values
    V1 = V(:,1:k);                 % dominant right singular vectors
    V2 = V(:,k+1:end);             % neglected right singular vectors
    Omega1 = V1'*Omega;
    Omega2 = V2'*Omega;
                                   % compute error bound (2-norm)
    dat(k,2)= sqrt(norm(Sig2)^2+norm(Sig2*Omega2*pinv(Omega1))^2);
 end

% plot error

 figure(1), clf
 semilogy([1:kmax],dat(:,1),'k.-','linewidth',2,'markersize',20)
 hold on
 semilogy([1:kmax],dat(:,2),'r.-','linewidth',2,'markersize',20)

 title(sprintf('H/M/T Theorem 9.1 with $m=%d$, $n=%d$', ...
    m,n),'interpreter','latex','fontsize',18)
 leg = legend('true error', 'Theorem 9.1','location','northeast');
 set(leg,'fontsize',16,'interpreter','latex')
 xlabel('$k$','fontsize',18,'interpreter','latex')
