function [W, obj] = lsd(A, r, max_iter, check_step, verbose, W0, nSeed)
% Left-Stochastic Matrix Decomposition (LSD)
% He Zhang, Sep. 2012
% Ref: R. Arora et al.:Clustering by left-stochastic matrix factorization, ICML 2011 

n = size(A,1);

if ~exist('verbose', 'var') || isempty(verbose)
    verbose = true;
end

if ~exist('nSeed', 'var') || isempty(nSeed)
    nSeed = 0;
end

if ~exist('W0', 'var') || isempty(W0)
    rand('seed',nSeed);
    W0 = rand(n,r);
end
W0 = bsxfun(@rdivide, W0, sum(W0,2)+eps);

if n<8000
    c = abs(sum(sum(pinv(full(A)))))/r;
else
    [E,D,~] = svds(A,30);
    M = bsxfun(@rdivide, E, sqrt(diag(D))'+eps);
    c = norm(sum(M))^2/r;
end

W = W0;
if verbose
    obj = ComputeObj(c, A, W);
    fprintf('init obj=%.6f\n', obj);
end


for iter=1:max_iter
    if mod(iter,check_step)==0
        W_old = W;
        if verbose
            fprintf('iter=%d, ', iter);
        end
    end
    
    gradn = c*A*W;
    gradp = W*(W'*W);
    a = sum(W./(gradp+eps),2);
    b = sum(W.*gradn./(gradp+eps),2);
    W = W .* (bsxfun(@times, gradn, a) + 1) ./ (bsxfun(@plus, bsxfun(@times, gradp, a), b)+eps);
    
    if mod(iter,check_step)==0
        if verbose
            obj = ComputeObj(c, A, W);
            fprintf('diffW=%.12f, ', norm(W-W_old));
            fprintf('stoc_err=%.12f, ', norm(sum(W,2)-ones(n,1)));
            fprintf('obj=%.6f, ', obj);
            fprintf('\n');
        end
    end
    
    if mod(iter, check_step)==0
        obj = ComputeObj(c, A, W);
        fprintf('LSD: iter=%d, obj=%6.5f\n', iter, obj);
    end
    
end

if nargout>1
    obj = ComputeObj(c, A, W);
end

function obj = ComputeObj(c, A, W)
cA = c*A;
obj = norm(cA,'fro')^2 - 2 * trace(W'*cA*W) + norm(W'*W,'fro')^2;
