[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

unsupervised_decomposition.hxx
1/************************************************************************/
2/* */
3/* Copyright 2008-2011 by Michael Hanselmann and Ullrich Koethe */
4/* */
5/* This file is part of the VIGRA computer vision library. */
6/* The VIGRA Website is */
7/* http://hci.iwr.uni-heidelberg.de/vigra/ */
8/* Please direct questions, bug reports, and contributions to */
9/* ullrich.koethe@iwr.uni-heidelberg.de or */
10/* vigra@informatik.uni-hamburg.de */
11/* */
12/* Permission is hereby granted, free of charge, to any person */
13/* obtaining a copy of this software and associated documentation */
14/* files (the "Software"), to deal in the Software without */
15/* restriction, including without limitation the rights to use, */
16/* copy, modify, merge, publish, distribute, sublicense, and/or */
17/* sell copies of the Software, and to permit persons to whom the */
18/* Software is furnished to do so, subject to the following */
19/* conditions: */
20/* */
21/* The above copyright notice and this permission notice shall be */
22/* included in all copies or substantial portions of the */
23/* Software. */
24/* */
25/* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26/* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27/* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28/* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29/* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30/* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31/* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32/* OTHER DEALINGS IN THE SOFTWARE. */
33/* */
34/************************************************************************/
35
36
37#ifndef VIGRA_UNSUPERVISED_DECOMPOSITION_HXX
38#define VIGRA_UNSUPERVISED_DECOMPOSITION_HXX
39
40#include <numeric>
41#include "mathutil.hxx"
42#include "matrix.hxx"
43#include "singular_value_decomposition.hxx"
44#include "random.hxx"
45
46namespace vigra
47{
48
49/** \addtogroup Unsupervised_Decomposition Unsupervised Decomposition
50
51 Unsupervised matrix decomposition methods.
52**/
53//@{
54
55/*****************************************************************/
56/* */
57/* principal component analysis (PCA) */
58/* */
59/*****************************************************************/
60
61 /** \brief Decompose a matrix according to the PCA algorithm.
62
63 This function implements the PCA algorithm (principal component analysis).
64
65 \arg features must be a matrix with shape <tt>(numFeatures * numSamples)</tt>, which is
66 decomposed into the matrices
67 \arg fz with shape <tt>(numFeatures * numComponents)</tt> and
68 \arg zv with shape <tt>(numComponents * numSamples)</tt>
69
70 such that
71 \f[
72 \mathrm{features} \approx \mathrm{fz} * \mathrm{zv}
73 \f]
74 (this formula requires that the features have been centered around the mean by
75 <tt>\ref linalg::prepareRows&nbsp;(features, features, ZeroMean)</tt>).
76
77 The shape parameter <tt>numComponents</tt> determines the complexity of
78 the decomposition model and therefore the approximation quality (if
79 <tt>numComponents == numFeatures</tt>, the representation becomes exact).
80 Intuitively, <tt>fz</tt> is a projection matrix from the reduced space
81 into the original space, and <tt>zv</tt> is the reduced representation
82 of the data, using just <tt>numComponents</tt> features.
83
84 <b>Declaration:</b>
85
86 <b>\#include</b> <vigra/unsupervised_decomposition.hxx>
87
88 \code
89 namespace vigra {
90
91 template <class U, class C1, class C2, class C3>
92 void
93 principalComponents(MultiArrayView<2, U, C1> const & features,
94 MultiArrayView<2, U, C2> fz,
95 MultiArrayView<2, U, C3> zv);
96 }
97 \endcode
98
99 <b>Usage:</b>
100 \code
101 Matrix<double> data(numFeatures, numSamples);
102 ... // fill the input matrix
103
104 int numComponents = 3;
105 Matrix<double> fz(numFeatures, numComponents),
106 zv(numComponents, numSamples);
107
108 // center the data
109 prepareRows(data, data, ZeroMean);
110
111 // compute the reduced representation
112 principalComponents(data, fz, zv);
113
114 Matrix<double> model = fz*zv;
115 double meanSquaredError = squaredNorm(data - model) / numSamples;
116 \endcode
117 */
118template <class T, class C1, class C2, class C3>
119void
123{
124 using namespace linalg; // activate matrix multiplication and arithmetic functions
125
126 int numFeatures = rowCount(features);
127 int numSamples = columnCount(features);
128 int numComponents = columnCount(fz);
129 vigra_precondition(numSamples >= numFeatures,
130 "principalComponents(): The number of samples has to be larger than the number of features.");
131 vigra_precondition(numFeatures >= numComponents && numComponents >= 1,
132 "principalComponents(): The number of features has to be larger or equal to the number of components in which the feature matrix is decomposed.");
133 vigra_precondition(rowCount(fz) == numFeatures,
134 "principalComponents(): The output matrix fz has to be of dimension numFeatures*numComponents.");
135 vigra_precondition(columnCount(zv) == numSamples && rowCount(zv) == numComponents,
136 "principalComponents(): The output matrix zv has to be of dimension numComponents*numSamples.");
137
138 Matrix<T> U(numSamples, numFeatures), S(numFeatures, 1), V(numFeatures, numFeatures);
139 singularValueDecomposition(features.transpose(), U, S, V);
140
141 for(int k=0; k<numComponents; ++k)
142 {
143 rowVector(zv, k) = columnVector(U, k).transpose() * S(k, 0);
144 columnVector(fz, k) = columnVector(V, k);
145 }
146}
147
148/*****************************************************************/
149/* */
150/* probabilistic latent semantic analysis (pLSA) */
151/* see T Hofmann, Probabilistic Latent Semantic */
152/* Indexing for details */
153/* */
154/*****************************************************************/
155
156 /** \brief Option object for the \ref pLSA algorithm.
157 */
159{
160 public:
161 /** Initialize all options with default values.
162 */
164 : min_rel_gain(1e-4),
165 max_iterations(50),
166 normalized_component_weights(true)
167 {}
168
169 /** Maximum number of iterations which is performed by the pLSA algorithm.
170
171 default: 50
172 */
174 {
175 vigra_precondition(n >= 1,
176 "PLSAOptions::maximumNumberOfIterations(): number must be a positive integer.");
177 max_iterations = n;
178 return *this;
179 }
180
181 /** Minimum relative gain which is required for the algorithm to continue the iterations.
182
183 default: 1e-4
184 */
186 {
187 vigra_precondition(g >= 0.0,
188 "PLSAOptions::minimumRelativeGain(): number must be positive or zero.");
189 min_rel_gain = g;
190 return *this;
191 }
192
193 /** Normalize the entries of the zv result array.
194
195 If true, the columns of zv sum to one. Otherwise, they have the same
196 column sum as the original feature matrix.
197
198 default: true
199 */
201 {
202 normalized_component_weights = v;
203 return *this;
204 }
205
206 double min_rel_gain;
207 int max_iterations;
208 bool normalized_component_weights;
209};
210
211 /** \brief Decompose a matrix according to the pLSA algorithm.
212
213 This function implements the pLSA algorithm (probabilistic latent semantic analysis)
214 proposed in
215
216 T. Hofmann: <a href="http://www.cs.brown.edu/people/th/papers/Hofmann-UAI99.pdf">
217 <i>"Probabilistic Latent Semantic Analysis"</i></a>,
218 in: UAI'99, Proc. 15th Conf. on Uncertainty in Artificial Intelligence,
219 pp. 289-296, Morgan Kaufmann, 1999
220
221 \arg features must be a matrix with shape <tt>(numFeatures * numSamples)</tt> and
222 non-negative entries, which is decomposed into the matrices
223 \arg fz with shape <tt>(numFeatures * numComponents)</tt> and
224 \arg zv with shape <tt>(numComponents * numSamples)</tt>
225
226 such that
227 \f[
228 \mathrm{features} \approx \mathrm{fz} * \mathrm{zv}
229 \f]
230 (this formula applies when pLSA is called with
231 <tt>PLSAOptions.normalizedComponentWeights(false)</tt>. Otherwise, you must
232 normalize the features by calling <tt>\ref linalg::prepareColumns&nbsp;(features, features, UnitSum)</tt>
233 to make the formula hold).
234
235 The shape parameter <tt>numComponents</tt> determines the complexity of
236 the decomposition model and therefore the approximation quality.
237 Intuitively, features are a set of words, and the samples a set of
238 documents. The entries of the <tt>features</tt> matrix denote the relative
239 frequency of the words in each document. The components represents a
240 (presumably small) set of topics. The matrix <tt>fz</tt> encodes the
241 relative frequency of words in the different topics, and the matrix
242 <tt>zv</tt> encodes to what extend each topic explains the content of each
243 document.
244
245 The option object determines the iteration termination conditions and the output
246 normalization. In addition, you may pass a random number generator to pLSA()
247 which is used to create the initial solution.
248
249 <b>Declarations:</b>
250
251 <b>\#include</b> <vigra/unsupervised_decomposition.hxx>
252
253 \code
254 namespace vigra {
255
256 template <class U, class C1, class C2, class C3, class Random>
257 void
258 pLSA(MultiArrayView<2, U, C1> const & features,
259 MultiArrayView<2, U, C2> & fz,
260 MultiArrayView<2, U, C3> & zv,
261 Random const& random,
262 PLSAOptions const & options = PLSAOptions());
263
264 template <class U, class C1, class C2, class C3>
265 void
266 pLSA(MultiArrayView<2, U, C1> const & features,
267 MultiArrayView<2, U, C2> & fz,
268 MultiArrayView<2, U, C3> & zv,
269 PLSAOptions const & options = PLSAOptions());
270 }
271 \endcode
272
273 <b>Usage:</b>
274 \code
275 Matrix<double> words(numWords, numDocuments);
276 ... // fill the input matrix
277
278 int numTopics = 3;
279 Matrix<double> fz(numWords, numTopics),
280 zv(numTopics, numDocuments);
281
282 pLSA(words, fz, zv, PLSAOptions().normalizedComponentWeights(false));
283
284 Matrix<double> model = fz*zv;
285 double meanSquaredError = (words - model).squaredNorm() / numDocuments;
286 \endcode
287 */
288doxygen_overloaded_function(template <...> void pLSA)
289
290
291template <class U, class C1, class C2, class C3, class Random>
292void
293pLSA(MultiArrayView<2, U, C1> const & features,
296 Random const& random,
297 PLSAOptions const & options = PLSAOptions())
298{
299 using namespace linalg; // activate matrix multiplication and arithmetic functions
300
301 int numFeatures = rowCount(features);
302 int numSamples = columnCount(features);
303 int numComponents = columnCount(fz);
304 vigra_precondition(numFeatures >= numComponents && numComponents >= 1,
305 "pLSA(): The number of features has to be larger or equal to the number of components in which the feature matrix is decomposed.");
306 vigra_precondition(rowCount(fz) == numFeatures,
307 "pLSA(): The output matrix fz has to be of dimension numFeatures*numComponents.");
308 vigra_precondition(columnCount(zv) == numSamples && rowCount(zv) == numComponents,
309 "pLSA(): The output matrix zv has to be of dimension numComponents*numSamples.");
310
311 // random initialization of result matrices, subsequent normalization
312 UniformRandomFunctor<Random> randf(random);
313 initMultiArray(destMultiArrayRange(fz), randf);
314 initMultiArray(destMultiArrayRange(zv), randf);
315 prepareColumns(fz, fz, UnitSum);
316 prepareColumns(zv, zv, UnitSum);
317
318 // init vars
319 double eps = 1.0/NumericTraits<U>::max(); // epsilon > 0
320 double lastChange = NumericTraits<U>::max(); // infinity
321 double err = 0;
322 double err_old;
323 int iteration = 0;
324
325 // expectation maximization (EM) algorithm
326 Matrix<U> columnSums(1, numSamples);
327 features.sum(columnSums);
328 Matrix<U> expandedSums = ones<U>(numFeatures, 1) * columnSums;
329
330 while(iteration < options.max_iterations && (lastChange > options.min_rel_gain))
331 {
332 Matrix<U> fzv = fz*zv;
333
334 //if(iteration%25 == 0)
335 //{
336 //std::cout << "iteration: " << iteration << std::endl;
337 //std::cout << "last relative change: " << lastChange << std::endl;
338 //}
339
340 Matrix<U> factor = features / pointWise(fzv + (U)eps);
341 zv *= (fz.transpose() * factor);
342 fz *= (factor * zv.transpose());
343 prepareColumns(fz, fz, UnitSum);
344 prepareColumns(zv, zv, UnitSum);
345
346 // check relative change in least squares model fit
347 Matrix<U> model = expandedSums * pointWise(fzv);
348 err_old = err;
349 err = (features - model).squaredNorm();
350 //std::cout << "error: " << err << std::endl;
351 lastChange = abs((err-err_old) / (U)(err + eps));
352 //std::cout << "lastChange: " << lastChange << std::endl;
353
354 iteration += 1;
355 }
356 //std::cout << "Terminated after " << iteration << " iterations." << std::endl;
357 //std::cout << "Last relative change was " << lastChange << "." << std::endl;
358
359 if(!options.normalized_component_weights)
360 {
361 // undo the normalization
362 for(int k=0; k<numSamples; ++k)
363 columnVector(zv, k) *= columnSums(0, k);
364 }
365}
366
367template <class U, class C1, class C2, class C3>
368inline void
369pLSA(MultiArrayView<2, U, C1> const & features,
370 MultiArrayView<2, U, C2> & fz,
371 MultiArrayView<2, U, C3> & zv,
372 PLSAOptions const & options = PLSAOptions())
373{
374 RandomNumberGenerator<> generator(RandomSeed);
375 pLSA(features, fz, zv, generator, options);
376}
377
378//@}
379
380} // namespace vigra
381
382
383#endif // VIGRA_UNSUPERVISED_DECOMPOSITION_HXX
384
Base class for, and view to, MultiArray.
Definition multi_array.hxx:705
U sum() const
Definition multi_array.hxx:1805
MultiArrayView< N, T, StridedArrayTag > transpose() const
Definition multi_array.hxx:1569
Option object for the pLSA algorithm.
Definition unsupervised_decomposition.hxx:159
PLSAOptions & minimumRelativeGain(double g)
Definition unsupervised_decomposition.hxx:185
PLSAOptions & normalizedComponentWeights(bool v=true)
Definition unsupervised_decomposition.hxx:200
PLSAOptions()
Definition unsupervised_decomposition.hxx:163
PLSAOptions & maximumNumberOfIterations(unsigned int n)
Definition unsupervised_decomposition.hxx:173
Definition random.hxx:785
Definition matrix.hxx:125
MultiArrayView< 2, vluae_type, StridedArrayTag > transpose() const
Linear algebra functions.
Definition eigensystem.hxx:50
void principalComponents(MultiArrayView< 2, T, C1 > const &features, MultiArrayView< 2, T, C2 > fz, MultiArrayView< 2, T, C3 > zv)
Decompose a matrix according to the PCA algorithm.
Definition unsupervised_decomposition.hxx:120
MultiArrayIndex columnCount(const MultiArrayView< 2, T, C > &x)
Definition matrix.hxx:684
void initMultiArray(...)
Write a value to every element in a multi-dimensional array.
void pLSA(...)
Decompose a matrix according to the pLSA algorithm.
unsigned int singularValueDecomposition(MultiArrayView< 2, T, C1 > const &A, MultiArrayView< 2, T, C2 > &U, MultiArrayView< 2, T, C3 > &S, MultiArrayView< 2, T, C4 > &V)
Definition singular_value_decomposition.hxx:75
MultiArrayView< 2, T, C > rowVector(MultiArrayView< 2, T, C > const &m, MultiArrayIndex d)
Definition matrix.hxx:697
MultiArrayIndex rowCount(const MultiArrayView< 2, T, C > &x)
Definition matrix.hxx:671
FFTWComplex< R >::NormType abs(const FFTWComplex< R > &a)
absolute value (= magnitude)
Definition fftw3.hxx:1002
FFTWComplex< R >::SquaredNormType squaredNorm(const FFTWComplex< R > &a)
squared norm (= squared magnitude)
Definition fftw3.hxx:1044
void prepareColumns(...)
Standardize the columns of a matrix according to given DataPreparationGoals.
MultiArrayView< 2, T, C > columnVector(MultiArrayView< 2, T, C > const &m, MultiArrayIndex d)
Definition matrix.hxx:727

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.12.2 (Mon Apr 14 2025)