newton_raphson.m

Go to the documentation of this file.

function [xnew, resnorm, residual, exitflag, output] = newton_raphson(funptr, x0, params)
% function [vnew, resnorm, residual, exitflags, output] = newton_raphson(funptr, vold, params)
% Global projected Levenberg-Marquard method

if ~isfield(params, 'beta')
  params.beta = 0.1;
end
if ~isfield(params, 'gamma')
  params.gamma = 0.999995;
end
if ~isfield(params, 'sigma')
  params.sigma = 1e-4;
end
if ~isfield(params, 'TolLineSearch')
  params.TolLineSearch = 1e-10;
end
if ~isfield(params, 'TolFuncCount')
  params.TolFuncCount = inf;
end
if ~isfield(params, 'TolRes')
  params.TolRes = 1e-8;
end
if ~isfield(params, 'TolChange')
  params.TolChange = 1e-10;
end
if ~isfield(params, 'maxIter')
  params.maxIter = inf;
end

if ~isfield(params, 'debug')
  params.debug = false;
end

xnew = [];
xold = x0;
k = 0;
Fold = inf;
stepsizes = [];
foms = [];
resmaxes = [];
imaxes = [];

[F, J] = funptr(xold);
%J = computed_jacobian(J, xold, funptr);
[n, m] = size(J);
f  = sum(F.*F);
nF = sqrt(f);
funcI = 1;

while true;

  % function converged to zero
  [max_F, max_i] = max(abs(F));
  imaxes = [imaxes, max_i];
  resmaxes = [resmaxes, max_F];
  disp(['i = ', num2str(imaxes(end)), ' max = ', num2str(resmaxes(end))]);
  if resmaxes(end) < params.TolRes
    disp(resmaxes(end));
    exitflag = 1;
    if isempty(xnew)
      xnew = xold;
    end
    disp('Found solution.');
    disp(resmaxes(end));
%    pause;
    break;
  end

  % change in residual too small
  if max(abs(F-Fold)) < params.TolChange
    exitflag = 3;
    disp('Change in residual too small');
    break;
  end

  %one_rows = find(xold(1:400) > 1-100*eps);
%  if ~isempty(one_rows)
%      keyboard;
%  end
  %  TI = repmat(one_rows', m, 1);
  %  TI = TI(:);
  %  TJ = repmat((1:m)', 1, length(one_rows));
  %  TJ = TJ(:);
  %  T  = sparse(TI, TJ, ones(1, length(one_rows)*m), n, m);

  %J(sub2ind([n, m], one_rows, one_rows)) = 0;
%  J(:, one_rows) = 0; %= %J - (J & T);

  gf = J'*F;
  defect = (J) \ -F;
  xnew = params.Px(xold + defect);

  % first order optimalities
  foms = [foms, max(gf)];

  [Fnew, Jnew] = funptr(xnew);
  fnew  = sum(Fnew.*Fnew);
  nFnew = sqrt(fnew);
  funcI = funcI + 1;

  if params.debug
    [OK, ok_mat] = check_jacobian(pxnext, funptr, Jcand);
  end

  if nFnew <= params.gamma * nF
    % found local iteration
    Fold = F;
    F    = Fnew;
    f    = fnew;
    nF   = nFnew;
    J    = Jnew;

  else
    tkc = params.beta;
    while true
      % line search too small
      if tkc < params.TolLineSearch
        exitflag = -4;
        tk = 0;
        disp('Line search failed.');
        break;
      end

      xnew  = params.Px(xold + tkc*defect);

      [Fnew, Jnew] = funptr(xnew);
      fnew  = sum(Fnew.*Fnew);
      nFnew = sqrt(fnew);
      funcI = funcI + 1;

      if sum(Fnew.*Fnew) <= f + params.sigma * gf' * defect
        Fold = F;
        F    = Fnew;
        f    = fnew;
        nF   = nFnew;

        tk   = tkc;
        break;
      else
        tkc = tkc * params.beta;
      end
    end
    if tk == 0
      break;
    end
  end

  stepsizes = [stepsizes, norm(xnew - xold)];
  % x change to small
  if max(abs(xnew - xold)) < params.TolChange
    exitflag = 2;
    disp('Step change too small.');
    break;
  end
  xold = xnew;

%  tempsols = [tempsols, xold];

  % too many iterations
  k    = k+1;
  if k > params.maxIter
    exitflag = 0;
    disp('Maximum number of iterations exceeded.');
    break;
  end

  if funcI > params.TolFuncCount
    exitflag = -1;
    disp('Maximum number of function evaluations exceeded.');
    break;
  end

end

resnorm              = f;
residual             = F;
output.iterations    = k;
output.funcCount     = funcI;
output.stepsize      = stepsizes;
output.firstorderopt = foms;
output.resmax        = resmaxes;

%if nargout == 6
%  tempsols = PCA_fixspace(tempsols, [], [], 2);
%end

end

function jac_comp = computed_jacobian(jac_test, X, fun)
  jac_comp = zeros(size(jac_test));

  for j = 1:size(jac_test, 2);
    addend = zeros(size(X));
    addend(j) = 1e-10;
    jac_comp(:,j) = 1/(2e-10) .* (fun(X + addend) - fun(X - addend));
  end
end
function [OK, ok_mat] = check_jacobian(Utest, fun, jac_test, epsilon)

  if nargin <= 3
    epsilon = 5e-2;
  end

  jac_comp = computed_jacobian(jac_test, Utest, fun);

  nz_elements = abs(jac_comp) > 1e6*eps;
  ok_mat = zeros(size(jac_test));
  ok_mat(nz_elements) = abs(jac_comp(nz_elements) - jac_test(nz_elements))./(abs(jac_comp(nz_elements))) > epsilon;
  ok_mat = ok_mat + (nz_elements) - (abs(full(jac_test)) > 100*eps);
  OK = all(~ok_mat(:));

  if ~OK
    disp('Jacobian incorrect: ');
    keyboard;
  end
end