cubic_search.m

Nonlinear dynamical factor analysis Matlab package

function [s, net, tnet, status] = cubic_search(...
    s0, net0, tnet0, basekls, restc, step, fs, tfs, data, ...
    params, t0, status, missing, notimedep),
% CUBIC_SEARCH  Line search by method of cubic approximation

% Copyright (C) 1999-2005 Antti Honkela, Harri Valpola,
% Xavier Giannakopoulos and Matti Tornio.
%
% This package comes with ABSOLUTELY NO WARRANTY; for details
% see License.txt in the program package.  This is free software,
% and you are welcome to redistribute it under certain conditions;
% see License.txt for details.

% Stopping criterion for the step size
epsilon = 1e-6;
% Stopping criterion for the cost function
epsilon2 = 1e-8;
% Maximum number of line search iterations
maxiters = 30;

fs = [];
tfs = [];

% Make sure the initial step size is large enough to avoid numerical problems
t0 = max(epsilon, t0);

% Extract model parameters
sdim = size(s0, 1);
N = size(s0, 2);

% Initial points for the line search
itercount0 = 3;
t1 = 0;
t2 = .5 * t0;
t3 = t0;

s1 = s0;
net1 = net0;
tnet1 = tnet0;

[s2, net2, tnet2] = update_s_and_net(s0, net0, tnet0, t2, step);
[s3, net3, tnet3] = update_s_and_net(s0, net0, tnet0, t3, step);

c1 = basekls;
c2 = restc + kldiv(tfs, fs, s2, data, net2, tnet2, params, ...
		   missing, notimedep, status, 1);
c3 = restc + kldiv(tfs, fs, s3, data, net3, tnet3, params, ...
		   missing, notimedep, status, 1);

t4 = nan;
c4 = nan;

if status.debug,
  fprintf('%.2f (%.4g) %.2f (%.4g) %.2f (%.4g)\n', c1, t1, c2, t2, c3, t3);
end

while ((c1 > c2) & (c2 > c3)),
  itercount0 = itercount0 + 1;

  tnew = (t1.^2 .* (c2-c3) + t2.^2 .* (c3-c1) + t3.^2 .* (c1-c2)) ./ ...
	 (2*(t1 .* (c2-c3) + t2 .* (c3-c1) + t3 .* (c1-c2)));
  
  %fprintf('New value for t3: old %f, new %f', 2*t3, 2*tnew);
  tnew = min([max([tnew .75*t3]) 5*t3]);
  %fprintf(' (%f)\n', 2*tnew);
  
  t4 = t2;
  c4 = c2;
  c2 = c3;
  s2 = s3;
  net2 = net3;
  tnet2 = tnet3;
  t2 = t3;
  %t3 = 2*t3;
  t3 = 2*tnew;
  
  [s3, net3, tnet3] = update_s_and_net(s0, net0, tnet0, t3, step);
  c3 = restc + kldiv(tfs, fs, s3, data, net3, tnet3, params, ...
		     missing, notimedep, status, 1);
  if status.debug,
    fprintf('%.2f (%.4g) %.2f (%.4g) %.2f (%.4g)\n', c1, t1, c2, t2, c3, t3);
  end
end

while (((c1 < c2) & (c2 < c3)) | ((c1 < c2) & (c2 > c3)) | ...
       ~isfinite(c2) | ~isfinite(c3)),
  itercount0 = itercount0 + 1;
  tnew = (t1.^2 .* (c2-c3) + t2.^2 .* (c3-c1) + t3.^2 .* (c1-c2)) ./ ...
	 (2*(t1 .* (c2-c3) + t2 .* (c3-c1) + t3 .* (c1-c2)));
  
  %fprintf('New value for t2: old %f, new %f', .5*t2, tnew);
  tnew = min([max([tnew .1*t2]) .5*t2]);
  %fprintf(' (%f)\n', tnew);
  
  t4 = t3;
  c4 = c3;
  c3 = c2;
  s3 = s2;
  net3 = net2;
  tnet3 = tnet2;
  t3 = t2;
  %t2 = .5*t2;
  t2 = tnew;
  
  [s2, net2, tnet2] = update_s_and_net(s0, net0, tnet0, t2, step);
  c2 = restc + kldiv(tfs, fs, s2, data, net2, tnet2, params, ...
		     missing, notimedep, status, 1);
  if status.debug,
    fprintf('%.2f (%.4g) %.2f (%.4g) %.2f (%.4g)\n', c1, t1, c2, t2, c3, t3);
  end
end

if ((c1 < c2) & (c2 > c3)),
  warning('Non-convex point configuration for line search');
  keyboard;
end

% Normal linesearch starts here
itercount = 0;
while (((t3 - t1) > epsilon) & ...
       ((abs(c2-c3) + abs(c3-c1) + abs(c1-c2)) > epsilon2)),

  tnew = (t1.^2 .* (c2-c3) + t2.^2 .* (c3-c1) + t3.^2 .* (c1-c2)) ./ ...
	 (2*(t1 .* (c2-c3) + t2 .* (c3-c1) + t3 .* (c1-c2)));

  % If more than 3 samples are available, use cubic interpolation
  if ~isnan(t4),
    tnew2 = cubic_minimum([t1 t2 t3 t4], [c1 c2 c3 c4]);
    % Check if the cubic approximation is valid, if not, revert
    % back to quadratic approximation
    if tnew2 > t1 & tnew2 < t3 & isreal(tnew2),
      tnew = tnew2;
    else
      warning('NDFA:LINESEARCH:cubinv', ...
	      'LINESEARCH: Cubic approximation invalid.');
    end
  end

  if (tnew == t2),
    warning('NDFA:LINESEARCH:neweqold', ...
	    'LINESEARCH: proposed new point equal to old midpoint.');
    if (t3 - t2) > (t2 - t1),
      tnew = t2 + .1 * (t3 - t2);
    else
      tnew = t2 - .1 * (t3 - t2);
    end
  end
  
  [snew, netnew, tnetnew] = update_s_and_net(s0, net0, tnet0, tnew, step);
  cnew = restc + kldiv(tfs, fs, snew, data, netnew, tnetnew, params, ...
		       missing, notimedep, status, 1);

  if (tnew > t2),
    if (cnew > c2),
      t4 = t3;
      c4 = c3;
      t3 = tnew;
      s3 = snew;
      net3 = netnew;
      tnet3 = netnew;
      c3 = cnew;
    else
      t4 = t1;
      c4 = c1;
      t1 = t2;
      s1 = s2;
      net1 = net2;
      tnet1 = tnet2;
      c1 = c2;
      t2 = tnew;
      s2 = snew;
      net2 = netnew;
      tnet2 = tnetnew;
      c2 = cnew;
    end
  else    % tnew < t2
    if (cnew > c2),
      t4 = t1;
      c4 = c1;
      t1 = tnew;
      s1 = snew;
      net1 = netnew;
      tnet1 = tnetnew;
      c1 = cnew;
    else
      t4 = t3;
      c4 = c3;
      t3 = t2;
      s3 = s2;
      net3 = net2;
      tnet3 = tnet2;
      c3 = c2;
      t2 = tnew;
      s2 = snew;
      net2 = netnew;
      tnet2 = tnetnew;
      c2 = cnew;
    end
  end

  if status.debug,
    fprintf('%.2f (%.4g) %.2f (%.4g) %.2f (%.4g)\n', c1, t1, c2, t2, c3, t3);
  end
  itercount = itercount + 1;

  if itercount > maxiters,
    fprintf('Line search failed, trouble...\n');
    break;
  end
end

Nh = size(status.kls, 2);
status.history.LS_iters(Nh) = itercount;
status.history.LS_iters0(Nh) = itercount0;

s = s2;
net = net2;
tnet = tnet2;
status.t0 = 2 * t2;

if status.debug,
  fprintf('Cost after mean update:      c=%f\n', c2);
end