function [f,D,E,Errs] = RR_torre_New(X,Y,eta,lambda,rho,gamma)
itrN = 10;
errX = zeros(itrN,1);
errD = zeros(itrN,1);
[d_x,n] = size(X);
d_y = size(Y,1);
I = eye(d_x+1);
D = X;
Dhat = [D;ones(1,n)];
T = Y * Dhat'*inv(Dhat*Dhat'+ gamma*eye(d_x+1));
E = X-D;
Z1 = X/norm(X);
Z2 = Dhat/norm(Dhat);
mu1 = d_x*n/4/norm(X(:),1);
mu2 = d_x*n/4/norm(D(:),1);
lambda = 1 / sqrt(min(d_x,n));
s = svd(X);
nunormX = sum(s);
eta = nunormX / norm(Y-T*Dhat,'fro');
for iter = 1:itrN
    disp([num2str(iter),'...']);
    T = Y * Dhat'*inv(Dhat*Dhat'+ gamma*eye(d_x+1));
    Dhat = (eta*T'*T+mu2*I)\(eta*T'*Y-Z2+mu2*[D;ones(1,n)]);
    beta = mu1 + mu2;
    Z = 1/beta * ( Z1 + mu1 * (X-E) + Z2(1:d_x,:) + mu2 * Dhat(1:d_x,:) );
    [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))';
    E = So( lambda / mu1 , X - D + Z1 / mu1 );
    Z1 = Z1 + mu1*(X-D-E);
    Z2 = Z2 + mu2*(Dhat-[D;ones(1,n)]);
    Xdif_ratio = norm(X-D-E,'fro')/norm(X,'fro');
    Ddif_ratio = norm(Dhat-[D;ones(1,n)],'fro')/norm(Dhat,'fro');
    errX(iter) = Xdif_ratio;
    errD(iter) = Ddif_ratio;
    if  Xdif_ratio <=1e-8 && Ddif_ratio<=1e-8
        break;
    end
    mu1 = min(rho*mu1,1e6);
    mu2 = min(rho*mu2,1e6);
    %eta = nunormX / norm(Y-T*Dhat,'fro')^2;
end
f = T;
Errs = {errX,errD};
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