function [f,D,E,Errs] = RR_torre(X,y,eta,lambda,rho,gamma)
itrN = 10;
errX = zeros(itrN,1);
errD = zeros(itrN,1);
[d,n] = size(X);
I = eye(d);
T = y*X'*(X*X'+gamma*I);
Dhat = pinv(T)*y;
D = Dhat;
E = X-D;
Z1 = X/norm(X);
Z2 = D/norm(D);
mu1 = d*n/norm(X(:),1);
mu2 = d*n/norm(D(:),1);
H = eye(size(D,2));
for iter = 1:itrN
    T = y*Dhat'/(Dhat*Dhat'+gamma*I);
    Dhat = (eta*T'*T+mu2*I)\(eta*T'*y-Z2+mu2*D*H);
    beta = norm(mu1*eye(n)+mu2*H*H',2);
    Z = D-1/beta*(-Z1+mu1*(D-(X-E))-Z2*H'+mu2*(D*H-Dhat)*H');
    [U,S,V] = svd(Z);
    s = diag(S);
    s = max(s-1/beta,0);
    D = U(:,1:length(s))*diag(s)*V(:,1:length(s))';
    Z = X-D+Z1/mu1;
    E = bsxfun(@times,max(sqrt(sum(Z.^2,1))-lambda/mu1,0),normc(Z));
    Z1 = Z1+mu1*(X-D-E);
    Z2 = Z2+mu2*(Dhat-D*H);
    Xdif_ratio = norm(X-D-E,'fro')/norm(X,'fro');
    Ddif_ratio = norm(Dhat-D*H,'fro')/norm(Dhat,'fro');
    errX(iter) = Xdif_ratio;
    errD(iter) = Ddif_ratio;
    if  Xdif_ratio<=1e-6 && Ddif_ratio<=1e-6
        break;
    end
    mu1 = min(rho*mu1,1e3);
    mu2 = min(rho*mu2,1e3);
end
f = T;
Errs = {errX,errD};