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