ndfa.m

Nonlinear dynamical factor analysis Matlab package

function [sources, net, tnet, params, status, missing, clamped, notimedep, fs, tfs] = ndfa(data, varargin),
% NDFA  Run nonlinear dynamical factor analysis
%
%  [sources, net, tnet, params, status, missing, clamped, notimedep, x_reco, s_pred] = NDFA(data, otherargs...)
%  runs the NDFA algorithm for given data (different samples in
%  different columns, different channels in different rows. Use
%  NaNs for missing values) and returns the resulting sources, 
%  net, tnet, parameters, status information and reconstructed 
%  observations.
%
%  If the number of return values expected is only one, the first
%  four values will be returned in a single stucture that can be
%  further passed as argument to the algorithm.
%
%  The additional arguments can be given as a structure with
%  the field name indicating the name of the argument, or as a
%  list as in NDFA(data, 'name1', val1, 'name2', val2, ...),
%  or as a combination of these as in
%  NDFA(data, struct, 'name1', val, 'name2', val2, ...).
%  In the latter case the values given in the list take precedence
%  over the ones specified in the structure.
%
% --------------------------------------------------------------------
%  ACCEPTED ARGUMENTS
%
%  The recognised arguments are as follows:
%
%  'sources'            The values of sources returned by previous a
%                       run of the algorithm.  This option can be used
%                       to continue a previous simulation.
%  'initsources'        Initial values for the sources to be used.
%  'searchsources'      Number of sources to use, to be initialised
%                       with PCA
%  'initclamped'        Initial values for clamped sources (use
%                       NaNs for missing values). This option can
%                       be used to supply control signals or
%                       otherwise known sources. The clamped
%                       sources are concatenated with the normal
%                       sources specified by 'sources',
%                       'initsources' or 'searchsources' argument.
%
%  'net'                The values of observation network weights
%                       returned by a previous run of the algorithm.
%                       This option can be used to continue a previous
%                       simulation.
%  'tnet'               The values of temporal network weights
%                       returned by a previous run of the algorithm.
%                       This option can be used to continue a previous
%                       simulation.
%  'hidneurons'         Number of hidden neurons to use in the
%                       observation MLP.
%  'thidneurons'        Number of hidden neurons to use in the
%                       temporal MLP.
%  'nonlin'             Nonlinear activation function to use in the
%                       MLPs.  Default value 'tanh' is currently the
%                       only supported value.
%
%  'params'             Structure of hyperparameter values returned by
%                       a previous run of the algorithm.  This option
%                       can be used to continue a previous simulation.
%
%  'status'             Structure of status information returned by a
%                       previous run of the algorithm.  This option
%                       can be used to continue a previous simulation.
%
%  'clamped'            The information on clamped sources returned by a
%                       previous run of the algorithm.  This option can be 
%                       used to continue a previous simulation.
%  'notimedep'          Sources which should have no time dependence with
%                       their successor. Can be used to partition the data
%                       into multiple parts with no time
%                       dependencies between.
%
%  'iters'              Number of iterations to run the algorithm.
%                       The default value is 300.
%                       Use value 0 to only evaluate the reconstructed
%                       observations and cost function value.
%  'embed'              Number of iterations during which embedded 
%                       sources are used.
%  'freeinitial'        Should the sources at first time instant have
%                       a zero mean prior (value 0) or a hierarchical
%                       prior automatically adapted to the correct
%                       value (value 1, default).
%  'approximation'      Nonlinearity approximation scheme to use.
%                       Supported values are 'hermite' (default) and
%                       'taylor'.
%  'updatealg'          Update algorithm.  Supported values are
%                       'old', 'grad', 'conjgrad'  and 
%                       'natgrad' (default).
%  'verbose'            Print out a summary of all simulation
%                       parameters at the beginning of simulation.
%                       The default value is 1 (yes), other possible
%                       values are 0 (no) and 2 (always).  With 1
%                       (yes).  The value is reset in the returned
%                       structure so further calls will not reprint
%                       the summary.
%  'debug'              Print out additional debug information for
%                       each iteration. Possible values are 0 (no
%                       debug information), 1 (normal information) 
%                       and 2 (extended information).
%  'epsilon'            Epsilon for stopping criterion.  The default
%                       value is 1e-6.
%  'nolearning'         Prevents updates to the networks or the 
%                       parameters, only sources are updated using 
%                       the preexisting mappings. Can be used to
%                       infer the sources with given observations
%                       and mappings.
%
% --------------------------------------------------------------------
%  RELATIONS BETWEEN ARGUMENTS
%
%  Exactly one of the arguments 'sources', 'initsources' and
%  'searchsources' must be set.
%
%  Either 'net' or 'hidneurons' must be set.  If 'net' is set,
%  'hidneurons' is not honoured.  Similarly either 'tnet' or
%  'thidneurons' must be set.  If 'tnet' is set, 'thidneurons' is not
%  honoured.
%
%
% --------------------------------------------------------------------
%  STOPPING CRITERIA
%
%  The iteration is stopped if the designated number of iterations is
%  reached, or the total number of iterations is greater than 400 and
%  a) the cost function increases for 10 iterations in a row; or
%  b) the cost function decreases by less than epsilon (see above)
%  on each iteration for 200 iterations in a row.
%
%
% --------------------------------------------------------------------
%  EXAMPLES
%
%  [sources, net, tnet, params, status] = ...
%      NDFA(data, 'searchsources', 5, 'hidneurons', 30, 'thidneurons', 20);
%            Extract 5 nonlinear factors from data using an
%            observation MLP with 30 hidden neurons and temporal MLP
%            with 20 hidden neurons and zero source initialisation.
%  
%  result = NDFA(data, 'initsources', my_s, 'hidneurons', 30, 'thidneurons', 20, 'iters', 50);
%            Extract nonlinear factors from data using an observation
%            MLP with 30 hidden neurons and temporal MLP with 20
%            hidden neurons and custom initialisation given by my_s
%            using 50 iterations of the algorithm.
%
%  result = NDFA(data, result, 'iters', 500);
%            Continue the previous simulation for 500 more iterations.
%
%  result = NDFA(data, 'searchsources', 6, 'initclamped', my_u, 'hidneurons', 30, 'thidneurons', 20);
%            Extract 6 nonlinear factors from data with clamped sources
%            my_u using an observation MLP with 30 hidden neurons and
%            temporal MLP with 20 hidden neurons.
%
%  result = NDFA(data, 'searchsources', 4, 'hidneurons', 30, 'thidneurons', 20, 'notimedep', 500);
%            Extract 4 nonlinear factors from two part data with the 
%            first part containing 500 samples using an observation 
%            MLP with 30 hidden neurons and a temporal MLP with 20 
%            hidden neurons.

% Copyright (C) 2002-2005 Harri Valpola, Antti Honkela 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.

% Read the arguments
args = readargs(varargin);

% Status
if ~isfield(args, 'status'),
  [status.iters, args] = getargdef(args, 'iters', 300);
  status.prune.hidneurons = 0;
  status.prune.thidneurons = 0;
  status.prune.sources = 0;
  status.prune.iters = 0;
  status.history = {};
  [status.embed.iters, args] = getargdef(args, 'embed', 0);
  [status.approximation, args] = ...
      getargdef(args, 'approximation', 'hermite');
  [status.updatealg, args] = getargdef(args, 'updatealg', 'natgrad');
  [status.epsilon, args] = getargdef(args, 'epsilon', 1e-6);
  [status.verbose, args] = getargdef(args, 'verbose', 1);
  [status.debug, args] = getargdef(args, 'debug', 0);
  [status.freeinitial, args] = getargdef(args, 'freeinitial', 1);
  
  status.kls = [];
  status.cputime = [cputime];
  status.version = 1.0;
  
  % How many iterations to wait before starting to update these values
  status.updatenet = 0;
  status.updatetnet = 0;
  if strcmp(status.updatealg, 'old'),
    [status.updatesrcs, args] = getargdef(args, 'updatesrcs', -50);
    [status.updatesrcvars, args] = getargdef(args, 'updatesrcvars', -50);
    [status.updateparams, args] = getargdef(args, 'updateparams', -100);
  else
    [status.updatesrcs, args] = getargdef(args, 'updatesrcs', -20);
    [status.updatesrcvars, args] = getargdef(args, 'updatesrcvars', -20);
    [status.updateparams, args] = getargdef(args, 'updateparams', -50);
  end
  status.t0 = 1;
else
  status = args.status;
  args = rmfield(args, 'status');
  % Support for obsolote models
  if ~isfield(status, 'version') || ...
     ~isnumeric(status.version) || ...
     (status.version < 1.01),
    status = convert_obsolote(status);
  end
  
  [status.iters, args] = getargdef(args, 'iters', status.iters);
  [status.embed.iters, args] = ...
      getargdef(args, 'embed', status.embed.iters);
  [status.prune.sources, args] = ...
      getargdef(args, 'prune', status.prune.sources);
  [status.prune.hidneurons, args] = ...
      getargdef(args, 'prune', status.prune.hidneurons);
  [status.prune.thidneurons, args] = ...
      getargdef(args, 'prune', status.prune.thidneurons);
  [status.prune.iters, args] = ...
      getargdef(args, 'prune', status.prune.iters);
  [status.approximation, args] = ...
      getargdef(args, 'approximation', status.approximation);
  [status.updatealg, args] = ...
      getargdef(args, 'updatealg', status.updatealg);
  [status.epsilon, args] = getargdef(args, 'epsilon', status.epsilon);
  [status.verbose, args] = getargdef(args, 'verbose', status.verbose);
  [status.debug, args] = getargdef(args, 'debug', status.debug);
end

if status.verbose,
  if isempty(status.kls),
    fprintf('Starting an NDFA simulation with parameters:\n');
  else
    fprintf('Continuing an NDFA simulation with parameters:\n');
    fprintf('Number of iterations so far: %d\n', length(status.kls));
  end
  fprintf('Number of signals: %d\n', size(data, 1));
  fprintf('Number of samples: %d\n', size(data, 2));
  switch status.approximation,
   case 'taylor'
    fprintf('Using Taylor approximation for the nonlinearity.\n');
   case 'hermite'
    fprintf('Using Gauss-Hermite approximation for the nonlinearity.\n');
   otherwise
    fprintf('Using an unknown approximation for the nonlinearity.\n');
  end
  switch status.updatealg,
   case 'old'
    fprintf('Using the old heuristic NDFA update algorithm.\n');
   case 'conjgrad'
    fprintf('Using the conjugate gradient update algorithm.\n');
   case 'natgrad'
    fprintf('Using the natural conjugate gradient update algorithm.\n');
   case 'grad'
    fprintf('Using a gradient update algorithm with line searches.\n');
   otherwise
    fprintf('Using an unknown update algorithm.\n');
  end
end

% Initialise the sources
if ~isfield(args, 'sources'),
  if isfield(args, 'initsources'),
    if status.verbose,
      fprintf('Using %d pre-initialised sources.\n', ...
 	      size(args.initsources, 1));
    end
    initsources = args.initsources;
    sources = acprobdist_alpha(...
	initsources, .0001 * ones(size(initsources)));
    args = rmfield(args, 'initsources');
  else
    if ~isfield(args, 'searchsources')
      error('Either sources, initsources or searchsources must be set!')
    end
    
    if status.verbose,
      fprintf('Using %d sources initialised with PCA.\n', ...
	      args.searchsources);
    end
    [data_em, s_mean] = embed(data, [], args.searchsources);
    zs = ones(size(s_mean));
    sources = acprobdist_alpha(s_mean, .0001*zs);
    status.embed.iters = -abs(status.embed.iters);
    if status.embed.iters,
      if status.verbose,
        fprintf('Embedding data.\n');
      end
      status.embed.datadim = size(data, 1);
      status.embed.timedim = size(data, 2);
      data = data_em;
    end
    args = rmfield(args, 'searchsources');
  end
else
  if status.verbose,
    fprintf('Using %d previously used sources.\n', size(args.sources, 1));
  end
  sources = args.sources;
  if status.embed.iters < 0,