function [ h ] = ssd( D,z,epsilon )
    [m, n]=size(D);
    D0=D; z0=z;
    g=pinv(D)*z;
    support=[];
    remain=1:n;
    while norm(z)>epsilon
        [~,l]=max(abs(g));
        support=[support,remain(l)];
        a=D(:,l);
        a_unit=a/(a'*a);
        z=z-(z'*a)*a_unit ;
        D=D-repmat(a'*D,m,1).*repmat(a_unit,1,size(D,2));
        D(:,l)=[];
        remain(l)=[];
        g=pinv(D)*z;
    end
    h=zeros(n,1);
    h(support)=pinv(D0(:,support))*z0;
end