% First-Order Potential reduction algorithm for solving the problem
%
%      minimize      (1/2)\|s+A'y-c\|^2
%      subject to     s >= 0
%
%      Input: A, c;
%      
%      Output: (y,s)
% 
%      V2: using a different stepsize for each coordinate
% 
[m,n]=size(A);
ee=ones(n,1);
toler=1.e-3;
rho=n+sqrt(n);
beta=0.5;
% 
AAT=inv(A*A');
%
% set initial solutions
%
iter=0;
%
% scale c vector to reduce large entry variables
%
c=(2/(1+max(abs(c))))*c;
%s=1+max(0,c);
s=max(1,sum(c)/n)*ee;
%
cs=c-s;
y=AAT*(A*cs);
r=A'*y-cs;
obj=norm(r)^2/2;
%
while (sqrt(2*obj)/(1+norm(c)) >= toler) & (iter < 20000),
   %
   % Compute the gradient vector
   % 
   g=(rho/obj)*r;
   % 
   % Scale the gradient
   %
   sg = s.*g-1;
   %
   % Further scaling
   %
   ds = min(beta./max(toler,sg), obj./(obj+rho*s.^2)).*sg;
   % 
   % Decide a simple step size
   %
   %theta  = max(ds);
   %if (theta < 0),
   %    theta=1;
   %end;
   %ads=A*ds;
   %aa=(2+rho/obj)*(ds'*ds)-(rho/obj)*(ads'*(AAT*ads));
   %alpha = sg'*ds/aa;
   %
   % Update s
   %
   s = s.*(1-ds);
   %
   % Update y
   %
   cs=c-s;
   y=AAT*(A*cs);
   r=A'*y-cs;
   %
   % Check the residual objective
   %
   obj=norm(r)^2/2;
   iter=iter+1;
end;
iter
norm(A'*y+s-c)/(1+norm(c))
%