function [V] = estimateS(X,num_train)

    XTX = X'*X;
    max_iter = 25;
    convg = zeros(1,max_iter);
    V = zeros(num_train, num_train);
    E = zeros(num_train, num_train);
    Y1 = zeros(num_train, num_train);
    Y2 = zeros(num_train, num_train);
    rho1 = 1;
    rho2 = 1;
    cc = 2;
    lambda = 1;
    mu0 = 0.1;
    for iter = 1 : max_iter
        Sraw = V + Y1/rho1;
        Shat = shrinkage(Sraw, mu0 / rho1);

        Vhat = (lambda*XTX + rho1*eye(num_train) + rho2*eye(num_train))\...
            (lambda*XTX + rho1*Shat -Y1 + rho2*E -Y2);

        Ehat = Vhat + Y2/rho2;
        Ehat(1 : num_train+1 : end)=0;

        E = Ehat;
        V = Vhat;
        S = Shat;

        Y1 = Y1 + rho1*(V-S);
        Y2 = Y2 + rho2*(V-E);

        rho1 = min(rho1*cc, 1e7);
        rho2 = min(rho2*cc, 1e7);

        Sdif_ratio = norm(S-V,'fro') / norm(S,'fro');
        Vdif_ratio = norm(V-E,'fro') / norm(V,'fro');
        convg(1,iter) = Sdif_ratio;
        if Sdif_ratio<1e-4 && Vdif_ratio<1e-4
            break;
        end
    end
end

