include/dai/emalg.h

   1 /*  This file is part of libDAI - http://www.libdai.org/
   2  *
   3  *  libDAI is licensed under the terms of the GNU General Public License version
   4  *  2, or (at your option) any later version. libDAI is distributed without any
   5  *  warranty. See the file COPYING for more details.
   6  *
   7  *  Copyright (C) 2009  Charles Vaske  [cvaske at soe dot ucsc dot edu]
   8  *  Copyright (C) 2009  University of California, Santa Cruz
   9  */
  10
  11
  12 #ifndef __defined_libdai_emalg_h
  13 #define __defined_libdai_emalg_h
  14
  15
  16 #include <vector>
  17 #include <map>
  18
  19 #include <dai/factor.h>
  20 #include <dai/daialg.h>
  21 #include <dai/evidence.h>
  22 #include <dai/index.h>
  23 #include <dai/properties.h>
  24
  25
  26 /// \file
  27 /// \brief Defines classes related to Expectation Maximization (EMAlg, ParameterEstimation, CondProbEstimation and SharedParameters)
  28 /// \todo Implement parameter estimation for undirected models / factor graphs.
  29
  30
  31 namespace dai {
  32
  33
  34 /// Base class for parameter estimation methods.
  35 /** This class defines the general interface of parameter estimation methods.
  36  *
  37  *  Implementations of this interface (see e.g. CondProbEstimation) should
  38  *  register a factory function (virtual constructor) via the static
  39  *  registerMethod() function.
  40  *  This factory function should return a pointer to a newly constructed
  41  *  object, whose type is a subclass of ParameterEstimation, and gets as
  42  *  input a PropertySet of parameters. After a subclass has been registered,
  43  *  instances of it can be constructed using the construct() method.
  44  *
  45  *  Implementations are responsible for collecting data from a probability
  46  *  vector passed to it from a SharedParameters container object.
  47  *
  48  *  The default registry only contains CondProbEstimation, named
  49  *  "CondProbEstimation".
  50  *
  51  *  \author Charles Vaske
  52  */
  53 class ParameterEstimation {
  54     public:
  55         /// Type of pointer to factory function.
  56         typedef ParameterEstimation* (*ParamEstFactory)( const PropertySet& );
  57
  58         /// Virtual destructor for deleting pointers to derived classes.
  59         virtual ~ParameterEstimation() {}
  60
  61         /// Virtual copy constructor.
  62         virtual ParameterEstimation* clone() const = 0;
  63
  64         /// General factory method that constructs the desired ParameterEstimation subclass
  65         /** \param method Name of the subclass that should be constructed;
  66          *  \param p Parameters passed to constructor of subclass.
  67          *  \note \a method should either be in the default registry or should be registered first using registerMethod().
  68          *  \throw UNKNOWN_PARAMETER_ESTIMATION_METHOD if the requested \a method is not registered.
  69          */
  70         static ParameterEstimation* construct( const std::string &method, const PropertySet &p );
  71
  72         /// Register a subclass so that it can be used with construct().
  73         static void registerMethod( const std::string &method, const ParamEstFactory &f ) {
  74             if( _registry == NULL )
  75                 loadDefaultRegistry();
  76             (*_registry)[method] = f;
  77         }
  78
  79         /// Estimate the factor using the accumulated sufficient statistics and reset.
  80         virtual Prob estimate() = 0;
  81
  82         /// Accumulate the sufficient statistics for \a p.
  83         virtual void addSufficientStatistics( const Prob &p ) = 0;
  84
  85         /// Returns the size of the Prob that should be passed to addSufficientStatistics.
  86         virtual size_t probSize() const = 0;
  87
  88     private:
  89         /// A static registry containing all methods registered so far.
  90         static std::map<std::string, ParamEstFactory> *_registry;
  91
  92         /// Registers default ParameterEstimation subclasses (currently, only CondProbEstimation).
  93         static void loadDefaultRegistry();
  94 };
  95
  96
  97 /// Estimates the parameters of a conditional probability table, using pseudocounts.
  98 /** \author Charles Vaske
  99  */
 100 class CondProbEstimation : private ParameterEstimation {
 101     private:
 102         /// Number of states of the variable of interest
 103         size_t _target_dim;
 104         /// Current pseudocounts
 105         Prob _stats;
 106         /// Initial pseudocounts
 107         Prob _initial_stats;
 108
 109     public:
 110         /// Constructor
 111         /** For a conditional probability \f$ P( X | Y ) \f$,
 112          *  \param target_dimension should equal \f$ | X | \f$
 113          *  \param pseudocounts are the initial pseudocounts, of length \f$ |X| \cdot |Y| \f$
 114          */
 115         CondProbEstimation( size_t target_dimension, const Prob &pseudocounts );
 116
 117         /// Virtual constructor, using a PropertySet.
 118         /** Some keys in the PropertySet are required.
 119          *  For a conditional probability \f$ P( X | Y ) \f$,
 120          *     - \a target_dimension should be equal to \f$ | X | \f$
 121          *     - \a total_dimension should be equal to \f$ |X| \cdot |Y| \f$
 122          *
 123          *  An optional key is:
 124          *     - \a pseudo_count which specifies the initial counts (defaults to 1)
 125          */
 126         static ParameterEstimation* factory( const PropertySet &p );
 127
 128         /// Virtual copy constructor
 129         virtual ParameterEstimation* clone() const { return new CondProbEstimation( _target_dim, _initial_stats ); }
 130
 131         /// Virtual destructor
 132         virtual ~CondProbEstimation() {}
 133
 134         /// Returns an estimate of the conditional probability distribution.
 135         /** The format of the resulting Prob keeps all the values for
 136          *  \f$ P(X | Y=y) \f$ in sequential order in the array.
 137          */
 138         virtual Prob estimate();
 139
 140         /// Accumulate sufficient statistics from the expectations in \a p
 141         virtual void addSufficientStatistics( const Prob &p );
 142
 143         /// Returns the required size for arguments to addSufficientStatistics().
 144         virtual size_t probSize() const { return _stats.size(); }
 145 };
 146
 147
 148 /// Represents a single factor or set of factors whose parameters should be estimated.
 149 /** To ensure that parameters can be shared between different factors during
 150  *  EM learning, each factor's values are reordered to match a desired variable
 151  *  ordering. The ordering of the variables in a factor may therefore differ
 152  *  from the canonical ordering used in libDAI. The SharedParameters
 153  *  class combines one or more factors (together with the specified orderings
 154  *  of the variables) with a ParameterEstimation object, taking care of the
 155  *  necessary permutations of the factor entries / parameters.
 156  *
 157  *  \author Charles Vaske
 158  */
 159 class SharedParameters {
 160     public:
 161         /// Convenience label for an index of a factor in a FactorGraph.
 162         typedef size_t FactorIndex;
 163         /// Convenience label for a grouping of factor orientations.
 164         typedef std::map<FactorIndex, std::vector<Var> > FactorOrientations;
 165
 166     private:
 167         /// Maps factor indices to the corresponding VarSets
 168         std::map<FactorIndex, VarSet> _varsets;
 169         /// Maps factor indices to the corresponding Permute objects that permute the canonical ordering into the desired ordering
 170         std::map<FactorIndex, Permute> _perms;
 171         /// Maps factor indices to the corresponding desired variable orderings
 172         FactorOrientations _varorders;
 173         /// Parameter estimation method to be used
 174         ParameterEstimation *_estimation;
 175         /// Indicates whether \c *this gets ownership of _estimation
 176         bool _ownEstimation;
 177
 178         /// Calculates the permutation that permutes the canonical ordering into the desired ordering
 179         /** \param varOrder Desired ordering of variables
 180          *  \param outVS Contains variables in \a varOrder represented as a VarSet
 181          *  \return Permute object for permuting variables in varOrder from the canonical libDAI ordering into the desired ordering
 182          */
 183         static Permute calculatePermutation( const std::vector<Var> &varOrder, VarSet &outVS );
 184
 185         /// Initializes _varsets and _perms from _varorders and checks whether their state spaces correspond with _estimation.probSize()
 186         void setPermsAndVarSetsFromVarOrders();
 187
 188     public:
 189         /// Constructor
 190         /** \param varorders  all the factor orientations for this parameter
 191          *  \param estimation a pointer to the parameter estimation method
 192          *  \param ownPE whether the constructed object gets ownership of \a estimation
 193          */
 194         SharedParameters( const FactorOrientations &varorders, ParameterEstimation *estimation, bool ownPE=false );
 195
 196         /// Construct a SharedParameters object from an input stream \a is and a factor graph \a fg
 197         /** \see \ref fileformats-emalg-sharedparameters
 198          *  \throw INVALID_EMALG_FILE if the input stream is not valid
 199          */
 200         SharedParameters( std::istream &is, const FactorGraph &fg );
 201
 202         /// Copy constructor
 203         SharedParameters( const SharedParameters &sp ) : _varsets(sp._varsets), _perms(sp._perms), _varorders(sp._varorders), _estimation(sp._estimation), _ownEstimation(sp._ownEstimation) {
 204             // If sp owns its _estimation object, we should clone it instead of copying the pointer
 205             if( _ownEstimation )
 206                 _estimation = _estimation->clone();
 207         }
 208
 209         /// Destructor
 210         ~SharedParameters() {
 211             // If we own the _estimation object, we should delete it now
 212             if( _ownEstimation )
 213                 delete _estimation;
 214         }
 215
 216         /// Collect the sufficient statistics from expected values (beliefs) according to \a alg
 217         /** For each of the relevant factors (that shares the parameters we are interested in),
 218          *  the corresponding belief according to \a alg is obtained and its entries are permuted
 219          *  such that their ordering corresponds with the shared parameters that we are estimating.
 220          *  Then, the parameter estimation subclass method addSufficientStatistics() is called with
 221          *  this vector of expected values of the parameters as input.
 222          */
 223         void collectSufficientStatistics( InfAlg &alg );
 224
 225         /// Estimate and set the shared parameters
 226         /** Based on the sufficient statistics collected so far, the shared parameters are estimated
 227          *  using the parameter estimation subclass method estimate(). Then, each of the relevant
 228          *  factors in \a fg (that shares the parameters we are interested in) is set according
 229          *  to those parameters (permuting the parameters accordingly).
 230          */
 231         void setParameters( FactorGraph &fg );
 232 };
 233
 234
 235 /// A MaximizationStep groups together several parameter estimation tasks (SharedParameters objects) into a single unit.
 236 /** \author Charles Vaske
 237  */
 238 class MaximizationStep {
 239     private:
 240         /// Vector of parameter estimation tasks of which this maximization step consists
 241         std::vector<SharedParameters> _params;
 242
 243     public:
 244         /// Default constructor
 245         MaximizationStep() : _params() {}
 246
 247         /// Construct MaximizationStep from a vector of parameter estimation tasks
 248         MaximizationStep( std::vector<SharedParameters> &maximizations ) : _params(maximizations) {}
 249
 250         /// Constructor from an input stream and a corresponding factor graph
 251         /** \see \ref fileformats-emalg-maximizationstep
 252          */
 253         MaximizationStep( std::istream &is, const FactorGraph &fg_varlookup );
 254
 255         /// Collect the beliefs from this InfAlg as expectations for the next Maximization step
 256         void addExpectations( InfAlg &alg );
 257
 258         /// Using all of the currently added expectations, make new factors with maximized parameters and set them in the FactorGraph.
 259         void maximize( FactorGraph &fg );
 260
 261     /// \name Iterator interface
 262     //@{
 263         /// Iterator over the parameter estimation tasks
 264         typedef std::vector<SharedParameters>::iterator iterator;
 265         /// Constant iterator over the parameter estimation tasks
 266         typedef std::vector<SharedParameters>::const_iterator const_iterator;
 267
 268         /// Returns iterator that points to the first parameter estimation task
 269         iterator begin() { return _params.begin(); }
 270         /// Returns constant iterator that points to the first parameter estimation task
 271         const_iterator begin() const { return _params.begin(); }
 272         /// Returns iterator that points beyond the last parameter estimation task
 273         iterator end() { return _params.end(); }
 274         /// Returns constant iterator that points beyond the last parameter estimation task
 275         const_iterator end() const { return _params.end(); }
 276     //@}
 277 };
 278
 279
 280 /// EMAlg performs Expectation Maximization to learn factor parameters.
 281 /** This requires specifying:
 282  *     - Evidence (instances of observations from the graphical model),
 283  *     - InfAlg for performing the E-step (which includes the factor graph),
 284  *     - a vector of MaximizationStep 's steps to be performed.
 285  *
 286  *  This implementation can perform incremental EM by using multiple
 287  *  MaximizationSteps. An expectation step is performed between execution
 288  *  of each MaximizationStep. A call to iterate() will cycle through all
 289  *  MaximizationStep 's. A call to run() will call iterate() until the
 290  *  termination criteria have been met.
 291  *
 292  *  Having multiple and separate maximization steps allows for maximizing some
 293  *  parameters, performing another E-step, and then maximizing separate
 294  *  parameters, which may result in faster convergence in some cases.
 295  *
 296  *  \author Charles Vaske
 297  */
 298 class EMAlg {
 299     private:
 300         /// All the data samples used during learning
 301         const Evidence &_evidence;
 302
 303         /// How to do the expectation step
 304         InfAlg &_estep;
 305
 306         /// The maximization steps to take
 307         std::vector<MaximizationStep> _msteps;
 308
 309         /// Number of iterations done
 310         size_t _iters;
 311
 312         /// History of likelihoods
 313         std::vector<Real> _lastLogZ;
 314
 315         /// Maximum number of iterations
 316         size_t _max_iters;
 317
 318         /// Convergence tolerance
 319         Real _log_z_tol;
 320
 321     public:
 322         /// Key for setting maximum iterations
 323         static const std::string MAX_ITERS_KEY;
 324         /// Default maximum iterations
 325         static const size_t MAX_ITERS_DEFAULT;
 326         /// Key for setting likelihood termination condition
 327         static const std::string LOG_Z_TOL_KEY;
 328         /// Default likelihood tolerance
 329         static const Real LOG_Z_TOL_DEFAULT;
 330
 331         /// Construct an EMAlg from several objects
 332         /** \param evidence Specifies the observed evidence
 333          *  \param estep Inference algorithm to be used for the E-step
 334          *  \param msteps Vector of maximization steps, each of which is a group of parameter estimation tasks
 335          *  \param termconditions Termination conditions @see setTermConditions()
 336          */
 337         EMAlg( const Evidence &evidence, InfAlg &estep, std::vector<MaximizationStep> &msteps, const PropertySet &termconditions )
 338           : _evidence(evidence), _estep(estep), _msteps(msteps), _iters(0), _lastLogZ(), _max_iters(MAX_ITERS_DEFAULT), _log_z_tol(LOG_Z_TOL_DEFAULT)
 339         {
 340               setTermConditions( termconditions );
 341         }
 342
 343         /// Construct an EMAlg from Evidence \a evidence, an InfAlg \a estep, and an input stream \a mstep_file
 344         /** \see \ref fileformats-emalg
 345          */
 346         EMAlg( const Evidence &evidence, InfAlg &estep, std::istream &mstep_file );
 347
 348         /// Change the conditions for termination
 349         /** There are two possible parameters in the PropertySet \a p:
 350          *    - \a max_iters maximum number of iterations
 351          *    - \a log_z_tol critical proportion of increase in logZ
 352          *
 353          *  \see hasSatisifiedTermConditions()
 354          */
 355         void setTermConditions( const PropertySet &p );
 356
 357         /// Determine if the termination conditions have been met.
 358         /** There are two sufficient termination conditions:
 359          *    -# the maximum number of iterations has been performed
 360          *    -# the ratio of logZ increase over previous logZ is less than the
 361          *       tolerance, i.e.,
 362          *       \f$ \frac{\log(Z_t) - \log(Z_{t-1})}{| \log(Z_{t-1}) | } < \mathrm{tol} \f$.
 363          */
 364         bool hasSatisfiedTermConditions() const;
 365
 366         /// Return the last calculated log likelihood
 367         Real logZ() const { return _lastLogZ.back(); }
 368
 369         /// Returns number of iterations done so far
 370         size_t Iterations() const { return _iters; }
 371
 372         /// Get the E-step method used
 373         const InfAlg& eStep() const { return _estep; }
 374
 375         /// Iterate once over all maximization steps
 376         /** \return Log-likelihood after iteration
 377          */
 378         Real iterate();
 379
 380         /// Iterate over a single MaximizationStep
 381         Real iterate( MaximizationStep &mstep );
 382
 383         /// Iterate until termination conditions are satisfied
 384         void run();
 385
 386     /// \name Iterator interface
 387     //@{
 388         /// Iterator over the maximization steps
 389         typedef std::vector<MaximizationStep>::iterator s_iterator;
 390         /// Constant iterator over the maximization steps
 391         typedef std::vector<MaximizationStep>::const_iterator const_s_iterator;
 392
 393         /// Returns iterator that points to the first maximization step
 394         s_iterator s_begin() { return _msteps.begin(); }
 395         /// Returns constant iterator that points to the first maximization step
 396         const_s_iterator s_begin() const { return _msteps.begin(); }
 397         /// Returns iterator that points beyond the last maximization step
 398         s_iterator s_end() { return _msteps.end(); }
 399         /// Returns constant iterator that points beyond the last maximization step
 400         const_s_iterator s_end() const { return _msteps.end(); }
 401     //@}
 402 };
 403
 404
 405 } // end of namespace dai
 406
 407
 408 /** \example example_sprinkler_em.cpp
 409  *  This example shows how to use the EMAlg class.
 410  */
 411
 412
 413 #endif