A small example of randmized singular value decompostion (SVD)

random_svd_example.m

clc;
close all;
clear all;

%% Creating a large dataset
measurement = 9000:3000:24000;
lp = length(measurement);
dtime = zeros(lp,1);
derror = zeros(lp,1);
rtime = zeros(lp,3);
rerror = zeros(lp,3);
for i=1:length(measurement)
    SS = randn(measurement(i),3000);
    tstart = tic;
    [U,E,V] = svd(SS);
    r = 0.1*size(SS,1); % target rank
    det_SVD = U(:,1:r)*E(1:r,1:r)*V(:,1:r)';
    dtime(i) = toc(tstart);
    derror(i) = 1/(size(SS,1)*size(SS,2))*norm(det_SVD - SS,'fro');
    %% Preprocessing
    % Select degree of power iteration
    index = 1;
    for q=3:3:9
    tstart = tic;
    % Create a random projection
    omega = normrnd(0,1/sqrt(size(SS,2)),size(SS,2),r);
    % Taking the random projection
    T = SS*omega;
    % Power iteration
    for iter = 1:q
        T = SS*(SS'*T);
    end
    
    %% Random SVD
    
    % QR Decomposition
    [Q,R] = qr(T,0); % economy QR decomposition
    % Project SS to Q
    gamma = Q'*SS;
    % randomized SVD
    [UU,E,V] = svd(gamma,'econ');
    U = Q*UU;
    rtime(i,index) = toc(tstart);
    
    %Reconstruct reduced rank data
    rand_SVD = Q*gamma;
    rerror(i,index) = 1/(size(SS,1)*size(SS,2))*norm(SS - rand_SVD,'fro');
    index = index + 1;
    end
end
set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','tex');
figure(1)
plot(measurement, sort(dtime))
xlabel('Number of Measurements','interpreter','latex')
ylabel('Runtime (in seconds)','interpreter','latex')
hold on
plot(measurement, sort(rtime(:,1)))
hold on
plot(measurement, sort(rtime(:,2)))
hold on
plot(measurement, sort(rtime(:,3)))
figure(2)
plot(measurement, sort(derror))
xlabel('Number of Measurements','interpreter','latex')
ylabel('Nomralized Approximation Error','interpreter','latex')
hold on
plot(measurement, sort(rerror(:,1)))
hold on
plot(measurement, sort(rerror(:,2)))
hold on
plot(measurement, sort(rerror(:,3)))