% (From  Lieven Vandenberghe)
% Input needed:  A,b
% e.g. m=500;n=100;A=randn(m,n);b=randn(m,1);
x=ones(n,1);
for k=1:50,
   y=A*x-b;
   val = sum(log(exp(y)+exp(-y)));
   grad = A'*((exp(y)-exp(-y))./(exp(y)+exp(-y)));
   if (norm(grad) < 1e-5), break; end;
   hess= 4*A'*diag(1./(exp(y)+exp(-y)).^2)*A;
   v=-hess\grad;
   t=1;
   tlogic=( ...
       sum(...
        log( exp( A*(x+t*v )-b)+exp( -A*(x+t*v)+b ) ) ) >...
                   (val +0.01*t*grad'*v)   );
   while  tlogic
     t=0.5*t;       % backtrack in linesearch
   tlogic=( ...
       sum(...
        log( exp( A*(x+t*v )-b)+exp( -A*(x+t*v)+b ) ) ) >...
                   (val +0.01*t*grad'*v)   );
   end;
   gvalue=sum( log( exp( A*(x+t*v )-b)+exp( -A*(x+t*v)+b ) )  );
   x=x+t*v;
   format long
   disp(['iteration: ',num2str(k),...
                 ' steplength: ',num2str(t),'  g value: ',num2str(gvalue)]);
   format short
end