function x=proj_polytope(y,A,b,N)
%INPUT
%================
% y ................. n-length vector
% A ................. mxn matrix
% b ................. m-length vector
% N ................ number of iterations
%OUTPUT
%================
% x ................. x=P_S(y), where S = {x: Ax<=b}

d=size(A);
m=d(1);
n=d(2);
lam=zeros(m,1);

L=2*max(eig(A*A'));
g=A*y-b;
for k=1:N
    lam=lam+2/L*(-A*(A'*lam)+g);
    lam=max(lam,0);
end
x=y-A'*lam;
