quadratic_search.m

非线型因素分析matlab仿真程序包

function [s, net, t2] = quadratic_search(s0, net0, basekls, restc, smstep, netmstep, fs, data, params, t0, status)
% QUADRATIC_SEARCH  Line search by method of quadratic approximation

% Copyright (C) 1999-2004 Antti Honkela, Harri Valpola,
% and Xavier Giannakopoulos.
%
% 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.

epsilon = 1e-6;
epsilon2 = 1e-8;
maxiters = 30;

sdim = size(s0, 1);
N = size(s0, 2);
t1 = 0;
t2 = .5 * t0;
t3 = t0;

s1 = s0;
net1 = net0;

[s2, net2] = update_s_and_net(s0, net0, t2, smstep, netmstep);
[s3, net3] = update_s_and_net(s0, net0, t3, smstep, netmstep);

c1 = basekls;
fs_tmp = feedfw(s2, net2, status.approximation);
c2 = kl_batch(fs_tmp{4}, s2, data, params) + restc + ...
     kl_static_split(net2, params, restc);
fs_tmp = feedfw(s3, net3, status.approximation);
c3 = kl_batch(fs_tmp{4}, s3, data, params) + restc + ...
     kl_static_split(net3, params, restc);

while ((c1 > c2) & (c2 > c3)),
  c2 = c3;
  s2 = s3;
  net2 = net3;
  t2 = t3;
  t3 = 2*t3;
  [s3, net3] = update_s_and_net(s0, net0, t3, smstep, netmstep);

  fs_tmp = feedfw(s3, net3, status.approximation);
  c3 = kl_batch(fs_tmp{4}, s3, data, params) + restc + ...
       kl_static_split(net3, params, restc);
  %fprintf('%.2f (%.4g) %.2f (%.4g) %.2f (%.4g)\n', c1, t1, c2, t2, c3, t3);
end

while (((c1 < c2) & (c2 < c3)) | ((c1 < c2) & (c2 > c3)) | ~isfinite(c2) | ~isfinite(c3)),
  c3 = c2;
  s3 = s2;
  net3 = net2;
  t3 = t2;
  t2 = .5*t2;
  [s2, net2] = update_s_and_net(s0, net0, t2, smstep, netmstep);
  fs_tmp = feedfw(s2, net2, status.approximation);
  c2 = kl_batch(fs_tmp{4}, s2, data, params) + restc + ...
       kl_static_split(net2, params, restc);
  %fprintf('%.2f (%.4g) %.2f (%.4g) %.2f (%.4g)\n', c1, t1, c2, t2, c3, t3);
end

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

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 (tnew == t2),
    warning('NLFA:LINESEARCH:neweqold', 'LINESEARCH: proposed new point equal to old midpoint -> trouble...');
    if (t3 - t2) > (t2 - t1),
      tnew = t2 + .1 * (t3 - t2);
    else
      tnew = t2 - .1 * (t3 - t2);
    end
  end
  
  [snew, netnew] = update_s_and_net(s0, net0, tnew, smstep, netmstep);
  fs_tmp = feedfw(snew, netnew, status.approximation);
  cnew = kl_batch(fs_tmp{4}, snew, data, params) + restc + ...
	 kl_static_split(netnew, params, restc);
  
  if (tnew > t2),
    if (cnew > c2),
      t3 = tnew;
      s3 = snew;
      net3 = netnew;
      c3 = cnew;
    else
      t1 = t2;
      s1 = s2;
      net1 = net2;
      c1 = c2;
      t2 = tnew;
      s2 = snew;
      net2 = netnew;
      c2 = cnew;
    end
  else    % tnew < t2
    if (cnew > c2),
      t1 = tnew;
      s1 = snew;
      net1 = netnew;
      c1 = cnew;
    else
      t3 = t2;
      s3 = s2;
      net3 = net2;
      c3 = c2;
      t2 = tnew;
      s2 = snew;
      net2 = netnew;
      c2 = cnew;
    end
  end
  %fprintf('%.2f (%.2f) %.2f (%.2f) %.2f (%.2f)\n', c1, t1, c2, t2, c3, t3);

  itercount = itercount + 1;
  if itercount > maxiters,
    [s, net, t2] = goldsect_search(s0, net0, basekls, smstep, netmstep, ...
				   fs, data, params, t2, status);
    return;
  end
end

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

%        z = 2.0*(x1 * (y2-y3) + x2 * (y3-y1) + x3 * (y1-y2))
%        if z == 0:
%            xnew = x2
%        else:
%            xnew = (x1**2 * (y2-y3) + x2**2 * (y3-y1) +
%                    x3**2 * (y1-y2)) / z
%        return xnew