function [x,cluster,v] = kmedian(data,k,InitialCentroids)

l = size(data,1);

% initialize centroids
x = InitialCentroids;

% initialize clusters
cluster = zeros(l,1);
for i = 1 : l
    d = inf;
    for j = 1 : k
        if norm(data(i,:)-x(j,:),1) < d
            d = norm(data(i,:)-x(j,:),1);
            cluster(i) = j;
        end
    end
end

% compute the objective function value
vold = 0;
for i = 1 : l
    vold = vold + norm(data(i,:)-x(cluster(i),:),1);
end

while true

    % update centroids
    for j = 1 : k
        ind = find(cluster == j);
        if ~isempty(ind)
            x(j,:) = median(data(ind,:));
        end
    end

    % update clusters
    for i = 1 : l
        d = inf;
        for j = 1 : k
            if norm(data(i,:)-x(j,:),1) < d
                d = norm(data(i,:)-x(j,:),1);
                cluster(i) = j;
            end
        end
    end

    % update objective function
    v = 0;
    for i = 1 : l
        v = v + norm(data(i,:)-x(cluster(i),:),1);
    end

    % stopping criterion
    if vold - v < 1e-5
        break
    else
        vold = v;
    end

end

end