function Svr = VRPCA(X,lambda1,numepochs,winit)
% function Svr = VRPCA(X,lambda1,numepochs)
% A basic implementation of the VRPCA algorithm for solving
% argmax_{w:||w||=1} w'*(X*X'/n)*w
% Input: 
% - X: d*n data matrix
% - lambda1: Leading eigenvalue of (X*X'/n) (used to measure optimization error for benchmarking); 
% - numepochs: Required number of epochs
% - winit: Initialization point (d*1 vector). If none is supplied, a random point is chosen
% Output: 
% - Svr: vector of optimization errors after each epoch

[d,n] = size(X);

%%%%%%%%%% Parameters %%%%%%%%%%%%%
m = n; % epoch length
eta = 1/(mean(sum(X.^2,1))*sqrt(n)); % step size
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Svr = zeros(numepochs,1);

% initialize w
if nargin<4
    winit = randn(d,1); 
end;
w = winit/norm(winit);
wrun = w;

disp('Starting VR-PCA')
for i=1:numepochs
    u = (1/n)*X*(X'*w);
    for j=1:m
        s = randi(n);
        wrun = wrun+eta*(((X(:,s)'*(wrun-w)))*X(:,s)+u);
        wrun = wrun/norm(wrun);
    end;
    w = wrun;
    Svr(i) = 1-((1/n)*sum((X'*w).^2))/lambda1;
    fprintf('Finished epoch %d, error %f\n',i,Svr(i))
end