% This is an implemetation of the first-order potential reduction algorithm
% for
%     min    \|Ax-b\|^2,     s.t.  x\ge 0
%
% Use (AA^T) inverse and further scaling (individual step-sizes) 
% Input data
%
[m,n]=size(A);
%
%an=[];
%for j=1:n;
%    an=[an; 1/norm(A(:,j))];
%end
%A=A*sparse(diag(an));
%
ee=ones(n,1);
toler=1.e-3;
rho=2*n+sqrt(n);
beta=0.5;
% 
AAT=inv(A*A');
%
% set initial solutions
%
iter=0;
%
% scale b vector to reduce large entry variables
%
b=(2/(1+max(abs(b))))*b;
d=A'*(AAT*b);
%
%x=1+max(0,d);
x=max(1,sum(d)/n)*ee;
r=A*x-b;
objo=norm(A*x-b);
r=A'*(AAT*r);
obj=norm(r)^2/2;
%
while (objo/(1+norm(b)) >= toler) & (iter < 30000),
   %
   % form -b term of the quadratic equation where
   % a=rho, c=-obj
    % Compute the gradient vector
   % 
   g=(rho/obj)*r;
   % 
   % Scale the gradient
   %
   dx = x.*g-1;
   %
   % Further scaling
   %
   dx = min(beta./max(toler,dx), obj./(obj+rho*x.^2)).*dx;
   % 
   % Decide a simple step size
   %
   %theta  = max(ds);
   %if (theta < 0),
   %    theta=1;
   %end;
   %alpha = beta/theta;
   %
   % Update x
   %
   x = x.*(1-dx);
   %
   % Check the residual objective
   %
   r=A*x-b;
   objo=norm(r);
   r=A'*(AAT*r);
   obj=norm(r)^2/2;
   iter=iter+1;
end;
iter
norm(A*x-b)/(1+norm(b))
%

