goldsect_search.m

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

function [s, net, t0] = goldsect_search(s0, net0, basekls, smstep, netmstep, fs, data, params, t0, status)
% GOLDSECT_SEARCH  Line search by golden section method

% 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;
alpha = 0.6180339887498949;

sdim = size(s0, 1);
N = size(s0, 2);
t1 = 0;
t2 = (1 - alpha) * 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);
[netc, restc] = kl_static_split(net2, params);
c2 = kl_batch(fs_tmp{4}, s2, data, params) + netc + 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 = t3 / (1 - alpha);
  [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 = t2 * (1 - alpha);
  [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

%l = t1 + (1-alpha) * (t3-t1);
l = t2;
m = t1 + alpha     * (t3-t1);

sl = s2;
netl = net2;
cl = c2;
[sm, netm] = update_s_and_net(s0, net0, m, smstep, netmstep);
fs_tmp = feedfw(sm, netm, status.approximation);
cm = kl_batch(fs_tmp{4}, sm, data, params) + restc + ...
     kl_static_split(netm, params, restc);

%fprintf('%.2f (%.4g) %.2f (%.4g) %.2f (%.4g) %.2f (%.4g)\n', c1, t1, cl, l, cm, m, c3, t3);

while (((t3 - t1) > epsilon)),

  if cl > cm,
    t1 = l;
    c1 = cl;
    l = m;
    sl = sm;
    netl = netm;
    m = t1 + alpha     * (t3-t1);
    cl = cm;
    [sm, netm] = update_s_and_net(s0, net0, m, smstep, netmstep);
    fs_tmp = feedfw(sm, netm, status.approximation);
    cm = kl_batch(fs_tmp{4}, sm, data, params) + restc + ...
	 kl_static_split(netm, params, restc);
  else
    t3 = m;
    c3 = cm;
    m = l;
    sm = sl;
    netm = netl;
    l = t1 + (1-alpha) * (t3-t1);
    cm = cl;
    [sl, netl] = update_s_and_net(s0, net0, l, smstep, netmstep);
    fs_tmp = feedfw(sl, netl, status.approximation);
    cl = kl_batch(fs_tmp{4}, sl, data, params) + restc + ...
	 kl_static_split(netl, params, restc);
  end

  %fprintf('%.2f (%.4g) %.2f (%.4g) %.2f (%.4g) %.2f (%.4g)\n', c1, t1, cl, l, cm, m, c3, t3);
    
end

if cl < cm,
  t0 = l;
  s = sl;
  net = netl;
else
  t0 = m;
  s = sm;
  net = netm;
end

%fs_tmp = feedfw(s, net);
%c = kl_batch(fs_tmp{4}, s, data, params) + kl_static(net, params);
%fprintf('After search: %f\n', c);