Rewrote implementation of response propagation in MR
[libdai.git] / 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