Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
A small example of randmized singular value decompostion (SVD)
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)))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment