function [D,Z,E]= LRSR_D(X,A,lambda,rho)

[d_x,n] = size(A);
[d_x,nt] = size(X);
Z = zeros(n,nt);
for j = 1:nt
    DifAZj = A - repmat(X(:,j),[1,n]);
    sqdis = sum(DifAZj.^2);
    [~,mini] = min(sqdis);
    Z(mini,j) = 1;
end
%Z = zeros(n,nt);
I = eye(d_x+1);
J = zeros(n,nt);
D = A*Z;
E = X-D;
Y1 = X/norm(X);
Y3 = Z/norm(Z);
%Y3 = zeros(n,nt);
mu1 = d_x*n/4/norm(X(:),1);
mu3 = 100*d_x*n/4/norm(D(:),1);
mu3 = 100*mu1;
s = svd(A);
nunormA = sum(s);
lambda1 = 1 / sqrt(min(d_x,n));
lambda2 = 0.01 / sqrt(min(d_x,n));
for iter = 1:10
    disp([num2str(iter),'...']);
    
    TempE =  X - D + Y1 / mu1; 
    for cn=1:nt
        Ecol = TempE(:,cn);
        if norm(Ecol) < lambda1 / mu1
            Ecol(:) = 0;
        end
        E(:,cn) = Ecol;
    end
    E = So( lambda1 / mu1 , X - D + Y1 / mu1 );
    
    tau_Z1= vrho(A'*A) + 0.000001;
    beta_Z1 = (mu1 + mu3) * tau_Z1 / 1;
    beta_Z = norm(mu1*A'*A+mu3*eye(n),2);
    FZk = A'*( -Y1 + mu1 * ( D - (X-E) ) ) ...
          + mu3 * (Z - J) + Y3;
    [U,S,V] = svd(Z - FZk / beta_Z);
    s = diag(S);
    s = max(s-1/beta_Z,0);
    Z = U(:,1:length(s))*diag(s)*V(:,1:length(s))';
    D = A*Z;
    J = max(Z + Y3/mu3 - lambda2/mu3,0);
    
    Y1 = Y1 + mu1*(X-D-E);
    Y3 = Y3 + mu3*(Z-J);
    Xdif_ratio = norm(X-D-E,'fro')/norm(X,'fro');
    Zdif_ratio = norm(Z-J,'fro')/norm(Z,'fro');
    if  Xdif_ratio <=1e-8
        break;
    end
    mu1 = min(rho*mu1,1e6);
    mu3 = min(rho*mu3,1e6);
end

end
function r = So(tau, X)
    % shrinkage operator
    r = sign(X) .* max(abs(X) - tau, 0);
end

function r = Do(tau, X)
    % shrinkage operator for singular values
    [U, S, V] = svd(X, 'econ');
    r = U*So(tau, S)*V';
end