EM bugfix. Convenience methods in Factor, Permute, Properties, EM.
[libdai.git] / include / dai / emalg.h
1 /* Copyright (C) 2009 Charles Vaske [cvaske at soe dot ucsc dot edu]
2 University of California Santa Cruz
3
4 This file is part of libDAI.
5
6 libDAI is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2 of the License, or
9 (at your option) any later version.
10
11 libDAI is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with libDAI; if not, write to the Free Software
18 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
19 */
20
21
22 #ifndef __defined_libdai_emalg_h
23 #define __defined_libdai_emalg_h
24
25
26 #include <vector>
27 #include <map>
28
29 #include <dai/factor.h>
30 #include <dai/daialg.h>
31 #include <dai/evidence.h>
32 #include <dai/index.h>
33 #include <dai/properties.h>
34
35
36 /// \file
37 /// \brief Defines classes related to Expectation Maximization: EMAlg, ParameterEstimation, CondProbEstimation and SharedParameters
38 /// \todo Describe EM file format
39
40
41 namespace dai {
42
43
44 /// Interface for a parameter estimation method.
45 /** This parameter estimation interface is based on sufficient statistics.
46 * Implementations are responsible for collecting data from a probability
47 * vector passed to it from a SharedParameters container object.
48 *
49 * Implementations of this interface should register a factory function
50 * via the static ParameterEstimation::registerMethod function.
51 * The default registry only contains CondProbEstimation, named
52 * "ConditionalProbEstimation".
53 */
54 class ParameterEstimation {
55 public:
56 /// A pointer to a factory function.
57 typedef ParameterEstimation* (*ParamEstFactory)( const PropertySet& );
58
59 /// Virtual destructor for deleting pointers to derived classes.
60 virtual ~ParameterEstimation() {}
61 /// Virtual copy constructor.
62 virtual ParameterEstimation* clone() const = 0;
63
64 /// General factory method for construction of ParameterEstimation subclasses.
65 static ParameterEstimation* construct( const std::string &method, const PropertySet &p );
66
67 /// Register a subclass so that it can be used with ParameterEstimation::construct.
68 static void registerMethod( const std::string &method, const ParamEstFactory &f ) {
69 if( _registry == NULL )
70 loadDefaultRegistry();
71 (*_registry)[method] = f;
72 }
73
74 /// Estimate the factor using the accumulated sufficient statistics and reset.
75 virtual Prob estimate() = 0;
76
77 /// Accumulate the sufficient statistics for p.
78 virtual void addSufficientStatistics( const Prob &p ) = 0;
79
80 /// Returns the size of the Prob that should be passed to addSufficientStatistics.
81 virtual size_t probSize() const = 0;
82
83 private:
84 /// A static registry containing all methods registered so far.
85 static std::map<std::string, ParamEstFactory> *_registry;
86
87 /// Registers default ParameterEstimation subclasses (currently, only CondProbEstimation).
88 static void loadDefaultRegistry();
89 };
90
91
92 /// Estimates the parameters of a conditional probability table, using pseudocounts.
93 class CondProbEstimation : private ParameterEstimation {
94 private:
95 /// Number of states of the variable of interest
96 size_t _target_dim;
97 /// Current pseudocounts
98 Prob _stats;
99 /// Initial pseudocounts
100 Prob _initial_stats;
101
102 public:
103 /// Constructor
104 /** For a conditional probability \f$ P( X | Y ) \f$,
105 * \param target_dimension should equal \f$ | X | \f$
106 * \param pseudocounts has length \f$ |X| \cdot |Y| \f$
107 */
108 CondProbEstimation( size_t target_dimension, const Prob &pseudocounts );
109
110 /// Virtual constructor, using a PropertySet.
111 /** Some keys in the PropertySet are required.
112 * For a conditional probability \f$ P( X | Y ) \f$,
113 * - target_dimension should be equal to \f$ | X | \f$
114 * - total_dimension should be equal to \f$ |X| \cdot |Y| \f$
115 *
116 * An optional key is:
117 * - pseudo_count, which specifies the initial counts (defaults to 1)
118 */
119 static ParameterEstimation* factory( const PropertySet &p );
120
121 /// Virtual copy constructor
122 virtual ParameterEstimation* clone() const { return new CondProbEstimation( _target_dim, _initial_stats ); }
123
124 /// Virtual destructor
125 virtual ~CondProbEstimation() {}
126
127 /// Returns an estimate of the conditional probability distribution.
128 /** The format of the resulting Prob keeps all the values for
129 * \f$ P(X | Y=y) \f$ in sequential order in the array.
130 */
131 virtual Prob estimate();
132
133 /// Accumulate sufficient statistics from the expectations in p.
134 virtual void addSufficientStatistics( const Prob &p );
135
136 /// Returns the required size for arguments to addSufficientStatistics.
137 virtual size_t probSize() const { return _stats.size(); }
138 };
139
140
141 /// A single factor or set of factors whose parameters should be estimated.
142 /** To ensure that parameters can be shared between different factors during
143 * EM learning, each factor's values are reordered to match a desired variable
144 * ordering. The ordering of the variables in a factor may therefore differ
145 * from the canonical ordering used in libDAI. The SharedParameters
146 * class couples one or more factors (together with the specified orderings
147 * of the variables) with a ParameterEstimation object, taking care of the
148 * necessary permutations of the factor entries / parameters.
149 */
150 class SharedParameters {
151 public:
152 /// Convenience label for an index into a factor in a FactorGraph.
153 typedef size_t FactorIndex;
154 /// Convenience label for a grouping of factor orientations.
155 typedef std::map<FactorIndex, std::vector<Var> > FactorOrientations;
156
157 private:
158 /// Maps factor indices to the corresponding VarSets
159 std::map<FactorIndex, VarSet> _varsets;
160 /// Maps factor indices to the corresponding Permute objects that permute the desired ordering into the canonical ordering
161 std::map<FactorIndex, Permute> _perms;
162 /// Maps factor indices to the corresponding desired variable orderings
163 FactorOrientations _varorders;
164 /// Parameter estimation method to be used
165 ParameterEstimation *_estimation;
166 /// Indicates whether the object pointed to by _estimation should be deleted upon destruction
167 bool _deleteEstimation;
168
169 /// Calculates the permutation that permutes the variables in varorder into the canonical ordering
170 static Permute calculatePermutation( const std::vector<Var> &varorder, VarSet &outVS );
171
172 /// Initializes _varsets and _perms from _varorders
173 void setPermsAndVarSetsFromVarOrders();
174
175 public:
176 /// Copy constructor
177 SharedParameters( const SharedParameters &sp );
178
179 /// Constructor
180 /** \param varorders all the factor orientations for this parameter
181 * \param estimation a pointer to the parameter estimation method
182 */
183 SharedParameters( const FactorOrientations &varorders, ParameterEstimation *estimation, bool deleteParameterEstimationInDestructor=0);
184
185 /// Constructor for making an object from a stream and a factor graph
186 SharedParameters( std::istream &is, const FactorGraph &fg_varlookup );
187
188 /// Destructor
189 ~SharedParameters() {
190 if( _deleteEstimation )
191 delete _estimation;
192 }
193
194 /// Collect the necessary statistics from expected values
195 void collectSufficientStatistics( InfAlg &alg );
196
197 /// Estimate and set the shared parameters
198 void setParameters( FactorGraph &fg );
199
200 void collectParameters( const FactorGraph& fg, std::vector< Real >& outVals, std::vector< Var >& outVarOrder );
201 };
202
203
204 /// A MaximizationStep groups together several parameter estimation tasks into a single unit.
205 class MaximizationStep {
206 private:
207 std::vector<SharedParameters> _params;
208
209 public:
210 /// Default constructor
211 MaximizationStep() : _params() {}
212
213 /// Constructor from a vector of SharedParameters objects
214 MaximizationStep( std::vector<SharedParameters> &maximizations ) : _params(maximizations) {}
215
216 /// Constructor from an input stream and a corresponding factor graph
217 MaximizationStep( std::istream &is, const FactorGraph &fg_varlookup );
218
219 /// Collect the beliefs from this InfAlg as expectations for the next Maximization step.
220 void addExpectations( InfAlg &alg );
221
222 /// Using all of the currently added expectations, make new factors with maximized parameters and set them in the FactorGraph.
223 void maximize( FactorGraph &fg );
224
225 /// @name iterator interface
226 //@{
227 typedef std::vector< SharedParameters >::iterator iterator;
228 typedef std::vector< SharedParameters >::const_iterator const_iterator;
229 iterator begin() { return _params.begin(); }
230 const_iterator begin() const { return _params.begin(); }
231 iterator end() { return _params.end(); }
232 const_iterator end() const { return _params.end(); }
233 //@}
234 };
235
236
237 /// EMAlg performs Expectation Maximization to learn factor parameters.
238 /** This requires specifying:
239 * - Evidence (instances of observations from the graphical model),
240 * - InfAlg for performing the E-step, which includes the factor graph,
241 * - a vector of MaximizationSteps steps to be performed.
242 *
243 * This implementation can perform incremental EM by using multiple
244 * MaximizationSteps. An expectation step is performed between execution
245 * of each MaximizationStep. A call to iterate() will cycle through all
246 * MaximizationSteps.
247 *
248 * Having multiple and separate maximization steps allows for maximizing some
249 * parameters, performing another E step, and then maximizing separate
250 * parameters, which may result in faster convergence in some cases.
251 */
252 class EMAlg {
253 private:
254 /// All the data samples used during learning
255 const Evidence &_evidence;
256
257 /// How to do the expectation step
258 InfAlg &_estep;
259
260 /// The maximization steps to take
261 std::vector<MaximizationStep> _msteps;
262
263 /// Number of iterations done
264 size_t _iters;
265
266 /// History of likelihoods
267 std::vector<Real> _lastLogZ;
268
269 /// Maximum number of iterations
270 size_t _max_iters;
271
272 /// Convergence tolerance
273 Real _log_z_tol;
274
275 public:
276 /// Key for setting maximum iterations @see setTermConditions
277 static const std::string MAX_ITERS_KEY;
278 /// Default maximum iterations @see setTermConditions
279 static const size_t MAX_ITERS_DEFAULT;
280 /// Key for setting likelihood termination condition @see setTermConditions
281 static const std::string LOG_Z_TOL_KEY;
282 /// Default likelihood tolerance @see setTermConditions
283 static const Real LOG_Z_TOL_DEFAULT;
284
285 /// Construct an EMAlg from all these objects
286 EMAlg( const Evidence &evidence, InfAlg &estep, std::vector<MaximizationStep> &msteps, const PropertySet &termconditions )
287 : _evidence(evidence), _estep(estep), _msteps(msteps), _iters(0), _lastLogZ(), _max_iters(MAX_ITERS_DEFAULT), _log_z_tol(LOG_Z_TOL_DEFAULT)
288 {
289 setTermConditions( termconditions );
290 }
291
292 /// Construct an EMAlg from an Evidence object, an InfAlg object, and an input stream
293 EMAlg( const Evidence &evidence, InfAlg &estep, std::istream &mstep_file );
294
295 /// Change the conditions for termination
296 /** There are two possible parameters in the PropertySet
297 * - max_iters maximum number of iterations
298 * - log_z_tol proportion of increase in logZ
299 *
300 * \see hasSatisifiedTermConditions()
301 */
302 void setTermConditions( const PropertySet &p );
303
304 /// Determine if the termination conditions have been met.
305 /** There are two sufficient termination conditions:
306 * -# the maximum number of iterations has been performed
307 * -# the ratio of logZ increase over previous logZ is less than the
308 * tolerance, i.e.,
309 * \f$ \frac{\log(Z_t) - \log(Z_{t-1})}{| \log(Z_{t-1}) | } < \mathrm{tol} \f$.
310 */
311 bool hasSatisfiedTermConditions() const;
312
313 /// Return the last calculated log likelihood
314 Real getLogZ() const { return _lastLogZ.back(); }
315
316 /// Returns number of iterations done so far
317 size_t getCurrentIters() const { return _iters; }
318
319 /// Get the iteration method used
320 const InfAlg& eStep() const { return _estep; }
321
322 /// Perform an iteration over all maximization steps
323 Real iterate();
324
325 /// Perform an iteration over a single MaximizationStep
326 Real iterate( MaximizationStep &mstep );
327
328 /// Iterate until termination conditions are satisfied
329 void run();
330
331 /// @name iterator interface
332 //@{ !!!
333 typedef std::vector< MaximizationStep >::iterator s_iterator;
334 typedef std::vector< MaximizationStep >::const_iterator const_s_iterator;
335 s_iterator s_begin() { return _msteps.begin(); }
336 const_s_iterator s_begin() const { return _msteps.begin(); }
337 s_iterator s_end() { return _msteps.end(); }
338 const_s_iterator s_end() const { return _msteps.end(); }
339 //@}
340 };
341
342
343 } // end of namespace dai
344
345
346 #endif