forwardselect.m

人工神经网络的源码编程

function [subset, H, l, U, A, w, P] = forwardSelect(F, Y, options)
% [subset, H, l, U, A, w, P]  = forwardSelect(F, Y, options)
%
% Regularised orthogonal least squares algorithm.
% See "Regularisation in the Selection of Radial Basis
% Function Centres", 1995, Orr, M.J.L., Neural Computation,
% 7(3):606-623.
%
% Inputs
%
%   F        Design matrix of selectable centres  (p-by-M)
%   Y	     Output training data                 (p-by-n)
%   options  Control options                      (string)
%
% Output
%
%   subset   Indices of selected columns of F     (1-by-m)
%   H        Subset of F                          (p-by-m)
%   l        regularisation parameter             (real and non-negative)
%   U        U the upper triangular tarnsform     (m-by-m)
%   A        inv(H'*H + l * U' * U)               (m-by-m)
%   w        A * H' * Y                           (m-by-n)
%   P        I - H * A * H'                       (p-by-p)
%

% Preliminaries.
[ph,M] = size(F);
[py,n] = size(Y);
if py ~= ph
  error('forwardSelect: design and outputs have incompatible dimensions')
end
p = py;

% Default options
Verbose = 0;
Flops = 0;
Global = 0;
Term = 'g';
MaxAge = 2;
MaxReg = 0;
ReEst = 'n';
Threshold = 0.9;
lambda = 0.1;
OLS = 0;

% process options
if nargin > 2

  % initialise
  i = 1;
  [arg, i] = getNextArg(options, i);

  % scan through arguments
  while ~isempty(arg)

    if strcmp(arg, '-v')

      % verbose output required
      Verbose = 1;

    elseif strcmp(arg, '-V')

      % verbose output required with compute cost reporting
      Verbose = 1;
      Flops = 1;

    elseif strcmp(arg, '-t')

      % specify termination criterion: the method must be specified
      [arg, i] = getNextArg(options, i);
      method_given = 1;
      if strcmp(lower(arg), 'var')
        % use fraction of explained variance to terminate
        Term = 'e';
      elseif strcmp(lower(arg), 'uev')
        % use UEV (unbiased expected variance) to terminate
        Term = 'u';
      elseif strcmp(lower(arg), 'fpe')
        % use FPE (final prediction error) to terminate
        Term = 'f';
      elseif strcmp(lower(arg), 'gcv')
        % use GCV (generalised cross-validation) to terminate
        Term = 'g';
      elseif strcmp(lower(arg), 'bic')
        % use BIC (Bayesian information criterion) to terminate
        Term = 'b';
      elseif strcmp(lower(arg), 'loo')
        % use GCV (generalised cross-validation) to terminate
        Term = 'o';
      else
        % no, the method wasn't specified after all
        method_given = 0;
      end
      if method_given
        % read next arg
        [arg, ii] = getNextArg(options, i);
      else
        fprintf('forwardSelect: terminate with VAR, UEV, FPE, GCV, BIC or LOO\n')
        error('forwardSelect: bad or missing argument for -t option')
      end

      % is the value given?
      th = str2num(arg);
      if ~isempty(th)
        % a value, threshold or maximun age, is specified
        good_value = 1;
        if Term == 'e' & th > 0 & th < 1
          % the value is a good one for threshold methods
          Threshold = th;
        elseif sum(Term == 'ufgbo') == 1 & th >= 1
          % the value is a good one for maximum age methods
          MaxAge = round(th);
        else
          % actually, it's not a good value
          good_value = 0;
        end
        if good_value
          % get ready to advance to next arg
          i = ii;
        else
          fprintf('forwardSelect: acceptable termination values are:\n')
          fprintf('  VAR: 0 < value < 1\n')
          fprintf('  UEV, FPE, GCV, BIC or LOO: value >= 1\n')
          error('forwardSelect: bad value for -t option')
        end
      end

    elseif strcmp(arg, '-m')

      % maximum number of regressors allowed in subset - needs value
      [arg, i] = getNextArg(options, i);
      mr = str2num(arg);
      good_value = 1;
      if ~isempty(mr)
        MaxReg = round(str2num(arg));
        if MaxReg < 1
          good_value = 0;
        end
      else
        % the value argument is mandatory
        good_value = 0;
      end
      if ~good_value
        fprintf('forwardSelect: positive maximum size required\n')
        error('forwardSelect: bad or missing value for -m option')
      end

    elseif strcmp(arg, '-g')

      % global regularisation required
      Global = 1;

      % is the initial value of lambda given?
      [arg, ii] = getNextArg(options, i);
      ll = str2num(arg);
      if ~isempty(ll)
        if ll < 0
          fprintf('forwardSelect: regularisation parameter should be > 0\n')
          error('forwardSelect: bad value for -g option')
        else
          lambda = ll;
          i = ii;
        end
      end

    elseif strcmp(arg, '-r')

      % turn global regularisation on
      Global = 1;

      % is method specified
      [arg, ii] = getNextArg(options, i);
      method_given = 1;
      if strcmp(lower(arg), 'uev')
        % use fraction of explained variance to re-estimate
        ReEst = 'u';
      elseif strcmp(lower(arg), 'fpe')
        % use FPE (final prediction error) to re-estimate
        ReEst = 'f';
      elseif strcmp(lower(arg), 'gcv')
        % use GCV (generalised cross-validation) to re-estimate
        ReEst = 'g';
      elseif strcmp(lower(arg), 'bic')
        % use BIC (Bayesian information criterion) to re-estimate
        ReEst = 'b';
      else
        % no, the method wasn't specified, set default
        method_given = 0;
        ReEst = 'g';
      end
      if method_given
        % advance to next arg and read it
        i = ii;
        [arg, ii] = getNextArg(options, i);
      end

      % is a value given?
      ll = str2num(arg);
      if ~isempty(ll)
        % an initial value for lambda is specified
        if ll >= 0
          lambda = ll;
          i = ii;
        else
          fprintf('forwardSelect: regularisation parameter should be > 0\n')
          error('forwardSelect: bad value for -r option')
        end
      end

    elseif strcmp(lower(arg), 'ols')

      % turn orthogonal least squares on
      OLS = 1;

    else

      fprintf('%s\n', options)
      for k = 1:i-length(arg)-1 fprintf(' '); end
      for k = 1:length(arg) fprintf('^'); end
      fprintf('\n')
      error('forwardSelect: unrecognised option')

    end

    % get next argument
    [arg, i] = getNextArg(options, i);

  end

end

% if global regularisation isn't being used then
% set lambda to zero
if Global == 0
  lambda = 0;
end

% initialise
Continue = 1;
YY = traceProduct(Y', Y);
m = 0;
AgeOfMin = 0;
if Verbose
  fprintf('\nforwardSelect\n')
  fprintf('pass   add    ')
  if ReEst ~= 'n'
    fprintf('  lambda  ')
  end
  if Term == 'e'
    fprintf('   VAR    ')
  elseif Term == 'u'
    fprintf('   UEV    ')
  elseif Term == 'f'
    fprintf('   FPE    ')
  elseif Term == 'g'
    fprintf('   GCV    ')
  elseif Term == 'b'
    fprintf('   BIC    ')
  else
    fprintf('   LOO    ')
  end
  if Flops
    fprintf('  flops   ')
    flops(0)
  end
  fprintf('\n')
end

% search for the most significant regressors.
while Continue

  % increment number of regressors
  m = m + 1;

  if m == 1 % first regressor - initialise

    % get error change associated with each regressor
    err = rowSum((F' * Y).^2) ./ (lambda + diagProduct(F', F));

    % select the maximum change
    mer = max(err);
    tot = mer / YY;
    can = find(err == mer);
    choose = can(1);
    subset = choose;

    % verbose
    if Verbose
      fprintf('%4d %5d    ', m, choose);
    end

    if OLS

      % initialise Ho, the current orthogonalised design matrix
      % and Hn, which is the same as Ho but with normalised columns
      % and some other useful stuff
      fj = F(:,choose);
      fjTfj = fj' * fj;
      Ho(:,1) = fj;
      diagHoTHo = fjTfj;
      HoTY = fj' * Y;
      Hn(:,1) = fj / fjTfj;

      % initialise Fm ready for second iteration
      Fm = F - Hn(:,1) * (Ho(:,1)' * F);

      % initialise upper triangular matrix Um
      Um = 1;

    else

      % initialise Hm
      fj = F(:,choose);
      Hm = fj;

      % initialise Fm ready for second iteration
      Fm = F;
      Fm = F - (fj / (fj' * fj)) * (fj' * F);

    end

  else

    % get error change associated with each regressor
    numerator = rowSum((Fm' * Y).^2);
    denominator = lambda + diagProduct(Fm', Fm);
    denominator(subset) = ones(length(subset),1); % avoid division by zero
    err = numerator ./ denominator;

    % select the maximum change
    mer = max(err);
    tot = tot + mer / YY;
    can = find(err == mer);
    choose = can(1);
    subset = [subset choose];

    % verbose
    if Verbose
      fprintf('%4d %5d    ', m, choose);
    end

    if OLS

      % collect next columns for Ho and Hn and other useful stuff
      fj = Fm(:,choose);
      fjTfj = fj' * fj;
      Ho(:,m) = fj;
      diagHoTHo = [diagHoTHo; fjTfj];
      HoTY = [HoTY; fj' * Y];
      Hn(:,m) = fj / fjTfj;

      % recompute Fm ready for next iteration
      Fm = Fm - Hn(:,m) * (Ho(:,m)' * Fm);

      % update Um
      Um = [Um Hn(:,1:m-1)' * F(:,choose); zeros(1,m-1) 1];

    else

      % update H
      Hm = [Hm F(:,choose)];

      % recompute Fm ready for next iteration
      fj = Fm(:,choose);
      Fm = Fm - (fj / (fj' * fj)) * (fj' * Fm);

    end

  end

  % Re-estimate lambda if required
  if ReEst ~= 'n'

    if OLS

      YTHoTHoY = diagProduct(HoTY,HoTY');
      ldiagHoTHo  = lambda + diagHoTHo;
      l2diagHoTHo = lambda + ldiagHoTHo;