gplmm.m

Go to the documentation of this file.

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

if ~isfield(params, 'beta')
  params.beta = 0.1;
end
if ~isfield(params, 'sigma')
  params.sigma = 1e-4;
end
if ~isfield(params, 'gamma')
  params.gamma = 0.999995;
end
if ~isfield(params, 'mu')
  params.mu = 1;
end
if ~isfield(params, 'TolFun')
  params.TolFun = 1e-10;
end
if ~isfield(params, 'TolLineSearch')
  params.TolLineSearch = 1e-10;
end
if ~isfield(params, 'TolFuncCount')
  params.TolFuncCount = 1000;
end
if ~isfield(params, 'TolRes')
  params.TolRes = 1e-8;
end
if ~isfield(params, 'TolChange')
  params.TolChange = 1e-13;
end
if ~isfield(params, 'maxIter')
  params.maxIter = 300;
end
if ~isfield(params, 'rho')
  params.rho = 1e-8;
end
if ~isfield(params, 'armijo_p')
  params.armijo_p = 2.1;
end
if isfield(params, 'HessFun')
  compute_hess_matrix = 0;
else
  compute_hess_matrix = 1;
end

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


tempsols = [];
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);

  % gradient f
  gf = J'*F;

  setup.type = 'nofill';
  setup.milu = 'row';
  setup.droptol = 1e-4;
  if isempty(xnew)
    if compute_hess_matrix
      Hess = J' * J;
%      [L,U] = ilu(Hess, setup);
%      dkU = gmres(@(X) Hess * X, -gf, 3, 1e-10, min(2500, m), L, U);
%      dkU = cgs(Hess, -gf, [], [], L, U);
      dkU = Hess \ (-gf);
    else
%      [HessTemp, precond] = params.HessFun(params, J, [], []);
      HessFunWrap = @(X) params.HessFun(params, J, [], X);
      dkU = bicgstab(HessFunWrap, -gf, 1e-12, 600);%,
%      precond.L,precond.U);
%      dkU = cgs(HessFunWrap, -gf, [], 300);
%      dkU = gmres(HessFunWrap, -gf, [], 1e-12, min(2500, m));%, precond.L, precond.U);
    end
    if condest(Hess) > 1e-12
      xnew = params.Px(xold + dkU);
    else
      xnew = xold;
    end
    Fold = funptr(xnew);
  end
  muk = min(params.mu * (Fold'*Fold), 1) / (k+1);

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

  if max(abs(Fold)) > params.TolRes
    if compute_hess_matrix
      Hess = J' * J + muk .* speye(m,m);
%      [L,U] = ilu(Hess, setup);
%      dkU = bicgstab(Hess, -gf, 1e-10, 600, L, U);
%      dkU = gmres(Hess, -gf, 2, 1e-10, min(2500, m), L, U);
      dkU = Hess \ (-gf);
    else
  %    [HessTemp, precond] = params.HessFun(params, J, [], []);
      HessFunWrap = @(X) params.HessFun(params, J, [], X) + muk .* X;
      dkU = gmres(HessFunWrap, -gf, 2, 1e-12, min(2500, m));%, precond.L, precond.U);
    end
  end
  pxnext = params.Px(xold + dkU);

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

  funcI = funcI + 1;
  fcand  = sum(Fcand.*Fcand);
  nFcand = sqrt(fcand);

  armijo_sk = pxnext - xold;
  if nFcand <= params.gamma * nF
    % found local iteration
    Fold = F;
    F    = Fcand;
    f    = fcand;
    nF   = nFcand;
    J    = Jcand;
    xnew = pxnext;

  else
    if gf' * armijo_sk <= - params.rho * norm(armijo_sk)^params.armijo_p
      % type line search direction
      sdir = armijo_sk;
%      disp('armijo');
    else
      % projected gradient line search direction
      sdir = -gf;
%      disp('gradient');
    end

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

      xktkc  = params.Px(xold + tkc*sdir);

      [Fxktkc, J] = funptr(xktkc);
%      J = computed_jacobian(J, xktkc, funptr);
      funcI = funcI + 1;

      if sum(Fxktkc.*Fxktkc) <= f + params.sigma * gf' * (xktkc - xold);
        Fold = F;
        F    = Fxktkc;
        f    = sum(F.*F);
        nF   = sqrt(f);

        xnew = xktkc;
        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