function x=lasso_GDNM(A,b,mu)

    ATA=A'*A;
    gamma=0.5/eigs(ATA,1);

    function [value,zeros_index]=prox_of_function(gamma,y,mu)
        value=[];
        zeros_index=[];
        for i=1:length(y)
            if y(i)>mu*gamma
                value=[value;y(i)-mu*gamma];
            elseif y(i)<-mu*gamma
                value=[value;y(i)+mu*gamma];
            else 
                value=[value;0];
                zeros_index=[zeros_index;i];
            end
        end
    end   

    Q=inv(eye(size(ATA))-gamma*ATA);
    P=Q-eye(size(Q));
    c=-gamma*Q*transpose(A)*b;

    function value=value_of_function(A,b,mu,x)
        value=0.5*norm(A*x-b)^2+mu*norm(x,1);
    end
    function value_of_phi=phi(A,b,gamma,mu,y)
        [prox_y,zeros_index]=prox_of_function(gamma,y,mu);
        value_of_phi=0.5*transpose(P*y)*y+transpose(c)*y+gamma*mu*norm(prox_y,1)+0.5*norm(y-prox_y)^2;
    end
    function grad=finding_gradient(A,b,gamma,mu,y)
        [prox_y,zeros_index]=prox_of_function(gamma,y,mu);
        grad=Q*y-prox_y+c;
    end

    function dk=finding_dk(A,b,gamma,mu,y)
        grad=finding_gradient(A,b,gamma,mu,y);
        [prox_y,zeros_index]=prox_of_function(gamma,y,mu);
        X=P;
        for i=zeros_index
            X(i,:)=Q(i,:);
        end
        dk=-linsolve(X,grad);
    end
    function tau=finding_tk(A,b,gamma,mu,y)
        tau=1;
        dk=finding_dk(A,b,gamma,mu,y);
        grad=finding_gradient(A,b,gamma,mu,y);
        while phi(A,b,gamma,mu,y+tau*dk)-phi(A,b,gamma,mu,y)-0.1*tau*grad'*dk>0
            tau=tau/2;
        end
    end

   y=zeros(length(c),1);
   iter=0;
   while iter<100
       y=y+finding_tk(A,b,gamma,mu,y)*finding_dk(A,b,gamma,mu,y);
       iter=iter+1
       value=0.5*norm(A*(Q*y+c)-b)^2+mu*norm(Q*y+c,1)
   end

   x=Q*y+c;
end

A=[4/7,3,6;12/7,2,-3;6/7,-6,2;0,24,0];
b=[104/49;347/49;-649/49;0];
mu=1/3;
x=lasso_GDNM(A,b,mu)