function [X, fn] = nlssubprob(A,T,Xinit,tol,maxiter)

% H, grad: output solution and gradient
% iter: #iterations used
% V, W: constant matrices
% Hinit: initial solution
% tol: stopping tolerance
% maxiter: limit of iterations
epsilon = 1e-10;
X = Xinit; %H is X here
%WtV = W'*V;
%WtW = W'*W; 

alpha = 1; beta = 0.1; sigmma = 0.01;
gradinit = (sum(A,1)  - T./(A*X+1e-12)'*A)';
gradinit_norm = norm(gradinit);
for iter=1:maxiter,  
  
  %X = X.*((T./(A*X)'*A)./(sum(A,1)+t))';
  %f = log_likelihood_MAP(A, T, X, t);
  
  
  %projgrad = norm(grad(grad < 0 | X >0));
  %if projgrad < tol*gradinit_norm,
  %  break
  %end
  

  
  %X = X + 1e-12;
  grad = (sum(A,1)  - T./(A*X+1e-12)'*A)';
  % search step size 
  for inner_iter=1:20,
  
    Xn = max(X - alpha*grad, 0); d = Xn-X;
    f  = log_likelihood_alspgrad_2( A, T, X);
    fn = log_likelihood_alspgrad_2( A, T, Xn);
    suff_decr = fn - f - sigmma*grad'*d < 0;
    %gradd=sum(sum(grad.*d)); dQd = sum(sum((WtW*d).*d));
    %suff_decr = 0.99*gradd + 0.5*dQd < 0;
    if inner_iter==1,
      decr_alpha = ~suff_decr; Xp = X;
    end
    if decr_alpha, 
      if suff_decr,
        X = Xn; break;
      else
        alpha = alpha * beta;
      end
    else
      if ~suff_decr | Xp == Xn,
        X = Xp; break;
      else
        alpha = alpha/beta; Xp = Xn;
      end
    end
  end
  
  if abs(fn - f) < epsilon
    %f = new_f;
    return;
  end
  
end

