SVD.m

%% Iterative SVD 
%  Calculate the largetst singular value and it's corrsponding
% left-singular vector and right-singular vector
% Author: chouclee
% Date: 10/01/2014
function [L,w,R] = SVD(S)
% S=[6,8,5,9;
%     8,5,4,6;
%     9,8,7,7;
%     6,7,3,8;
%     2,4,5,6;
%     9,8,9,6;
%     7,6,8,6;
%     7,8,6,8];
m = size(S,1); % number of rows
n = size(S,2); % number of columns

epsilon = 1e-6; % when sum of squared error is less than epsilon, stop iteration 

%% Generate random vector whose L2-norm equals to 1.
% generate a random vector, then normalize it
L = rand(m, 1);
L = L / sqrt(sum(L.^2, 1));

%% Iterate until ||S-w*L*R_transpose||=1
% initialize empty vector
R = zeros(1, n);

% initialize covergent flag
isConvergent = 0;

% variable to store sum of squared error of previous iteration
prevDelta = 0;
while ~isConvergent
    unnormalizedR = S'* L;                      % estimate unormalized R using L
    w = sqrt(sum(unnormalizedR.^2, 1));         % calculate L2-norm of unnormalized R
    R = unnormalizedR / w;                      % update normalized R
    unnormalizedL = S * R;                      % estimate unnormalized L using R
    w = sqrt(sum(unnormalizedL.^2, 1));         % calculate L2-norm of unnormalized L
    L = unnormalizedL / w;                      % update normalized L
    estimS = w * L * R';                        % update estimated S
    delta = sum(sum((estimS - S).^2, 1), 2);    % calculate sum of squared error
    if  abs(delta - prevDelta) < epsilon
        isConvergent = 1;
    end
    prevDelta = delta;
end