% big PCA demo. Barnett 10/11/21
clear
%M=1e4; N=3e3; % works
%M=3e4; N=1e4;   % big example
M=4e4; N=2e4;  % bigger
a = 0.9;   % sing val decay rate
K = min(N,ceil(log(1e16)/log(1/a)))   %   num rank
rng(1);  % fix seed
[V,~] = qr(randn(N,K),0);   % random real orthog
[U,~] = qr(randn(M,K),0);   % "
sig1 = 161.5259516482;     % exact  (note large to teach rel err)
sig = sig1 * a.^(0:K-1);   % true sing vals, decay, starting at sig1
A = U*(diag(sig)*V');      % form big dense, block outer

%rng(3);
%rng(0);
rng(0);
x=randn(N,1); x = x/norm(x);  % power iteration  (seed must differ from above!)
tic; nit=30;
for i=1:nit, x = A'*(A*x); n=norm(x); sig1est(i)=sqrt(n); x=x/norm(x); end
toc
(sig1est-sig1)/sig1;   % rel err conv

figure(1); clf;
semilogy(1:nit-1,abs(sig1est(1:end-1)/sig1 - 1),'r+-'); hold on;
semilogy(1:nit-1,abs(sig1est(1:end-1)/sig1est(end) - 1),'b+-');
xlabel('iteration $k$','interpreter','latex');
ylabel('est rel err $\varepsilon$','interpreter','latex');
hold on; plot(1:nit-1,0.05*a.^(4*(1:nit-1)),'m--');  % off by factor of 2 - why?
text(10,0.5,'Convergence plot $\varepsilon(k)$','interpreter','latex');
text(12,1e-2,'slope gives rate $\alpha$ or $\log a$','interpreter','latex','rotation',-38);
axis tight;
set(gcf,'paperposition',[0 0 3.6 3.6]);
print -dpng bigpca.png

% rand svd Halko-Martinsson-Tropp
tic;
Y = A*randn(N,30);
[Q,~]=qr(Y,0);
B = Q'*A;
S=svd(B);
toc
k=5;   % to test
S(1:k) ./ sig(1:k)'  % ratios to true

if 0
  tic; Strue = svd(A); toc   % cheat. 2sec for 1e4 * 3e3..  10 min for bigger
  k=5; (Strue(1:k)-sig(1:k)')/sig1    % rel err should be emach
end  

