[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
vigra/unsupervised_decomposition.hxx | ![]() |
00001 /************************************************************************/ 00002 /* */ 00003 /* Copyright 2008-2011 by Michael Hanselmann and Ullrich Koethe */ 00004 /* */ 00005 /* This file is part of the VIGRA computer vision library. */ 00006 /* The VIGRA Website is */ 00007 /* http://hci.iwr.uni-heidelberg.de/vigra/ */ 00008 /* Please direct questions, bug reports, and contributions to */ 00009 /* ullrich.koethe@iwr.uni-heidelberg.de or */ 00010 /* vigra@informatik.uni-hamburg.de */ 00011 /* */ 00012 /* Permission is hereby granted, free of charge, to any person */ 00013 /* obtaining a copy of this software and associated documentation */ 00014 /* files (the "Software"), to deal in the Software without */ 00015 /* restriction, including without limitation the rights to use, */ 00016 /* copy, modify, merge, publish, distribute, sublicense, and/or */ 00017 /* sell copies of the Software, and to permit persons to whom the */ 00018 /* Software is furnished to do so, subject to the following */ 00019 /* conditions: */ 00020 /* */ 00021 /* The above copyright notice and this permission notice shall be */ 00022 /* included in all copies or substantial portions of the */ 00023 /* Software. */ 00024 /* */ 00025 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */ 00026 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */ 00027 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */ 00028 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */ 00029 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */ 00030 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */ 00031 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */ 00032 /* OTHER DEALINGS IN THE SOFTWARE. */ 00033 /* */ 00034 /************************************************************************/ 00035 00036 00037 #ifndef VIGRA_UNSUPERVISED_DECOMPOSITION_HXX 00038 #define VIGRA_UNSUPERVISED_DECOMPOSITION_HXX 00039 00040 #include <numeric> 00041 #include "mathutil.hxx" 00042 #include "matrix.hxx" 00043 #include "singular_value_decomposition.hxx" 00044 #include "random.hxx" 00045 00046 namespace vigra 00047 { 00048 00049 /** \addtogroup Unsupervised_Decomposition Unsupervised Decomposition 00050 00051 Unsupervised matrix decomposition methods. 00052 **/ 00053 //@{ 00054 00055 /*****************************************************************/ 00056 /* */ 00057 /* principle component analysis (PCA) */ 00058 /* */ 00059 /*****************************************************************/ 00060 00061 /** \brief Decompose a matrix according to the PCA algorithm. 00062 00063 This function implements the PCA algorithm (principle component analysis). 00064 00065 \arg features must be a matrix with shape <tt>(numFeatures * numSamples)</tt>, which is 00066 decomposed into the matrices 00067 \arg fz with shape <tt>(numFeatures * numComponents)</tt> and 00068 \arg zv with shape <tt>(numComponents * numSamples)</tt> 00069 00070 such that 00071 \f[ 00072 \mathrm{features} \approx \mathrm{fz} * \mathrm{zv} 00073 \f] 00074 (this formula requires that the features have been centered around the mean by 00075 <tt>\ref linalg::prepareRows (features, features, ZeroMean)</tt>). 00076 00077 The shape parameter <tt>numComponents</tt> determines the complexity of 00078 the decomposition model and therefore the approximation quality (if 00079 <tt>numComponents == numFeatures</tt>, the representation becomes exact). 00080 Intuitively, <tt>fz</tt> is a projection matrix from the reduced space 00081 into the original space, and <tt>zv</tt> is the reduced representation 00082 of the data, using just <tt>numComponents</tt> features. 00083 00084 <b>Declaration:</b> 00085 00086 <b>\#include</b> <vigra/unsupervised_decomposition.hxx> 00087 00088 \code 00089 namespace vigra { 00090 00091 template <class U, class C1, class C2, class C3> 00092 void 00093 principleComponents(MultiArrayView<2, U, C1> const & features, 00094 MultiArrayView<2, U, C2> fz, 00095 MultiArrayView<2, U, C3> zv); 00096 } 00097 \endcode 00098 00099 <b>Usage:</b> 00100 \code 00101 Matrix<double> data(numFeatures, numSamples); 00102 ... // fill the input matrix 00103 00104 int numComponents = 3; 00105 Matrix<double> fz(numFeatures, numComponents), 00106 zv(numComponents, numSamples); 00107 00108 // center the data 00109 prepareRows(data, data, ZeroMean); 00110 00111 // compute the reduced representation 00112 principleComponents(data, fz, zv); 00113 00114 Matrix<double> model = fz*zv; 00115 double meanSquaredError = squaredNorm(data - model) / numSamples; 00116 \endcode 00117 */ 00118 template <class T, class C1, class C2, class C3> 00119 void 00120 principleComponents(MultiArrayView<2, T, C1> const & features, 00121 MultiArrayView<2, T, C2> fz, 00122 MultiArrayView<2, T, C3> zv) 00123 { 00124 using namespace linalg; // activate matrix multiplication and arithmetic functions 00125 00126 int numFeatures = rowCount(features); 00127 int numSamples = columnCount(features); 00128 int numComponents = columnCount(fz); 00129 vigra_precondition(numSamples >= numFeatures, 00130 "principleComponents(): The number of samples has to be larger than the number of features."); 00131 vigra_precondition(numFeatures >= numComponents && numComponents >= 1, 00132 "principleComponents(): The number of features has to be larger or equal to the number of components in which the feature matrix is decomposed."); 00133 vigra_precondition(rowCount(fz) == numFeatures, 00134 "principleComponents(): The output matrix fz has to be of dimension numFeatures*numComponents."); 00135 vigra_precondition(columnCount(zv) == numSamples && rowCount(zv) == numComponents, 00136 "principleComponents(): The output matrix zv has to be of dimension numComponents*numSamples."); 00137 00138 Matrix<T> U(numSamples, numFeatures), S(numFeatures, 1), V(numFeatures, numFeatures); 00139 singularValueDecomposition(features.transpose(), U, S, V); 00140 00141 for(int k=0; k<numComponents; ++k) 00142 { 00143 rowVector(zv, k) = columnVector(U, k).transpose() * S(k, 0); 00144 columnVector(fz, k) = columnVector(V, k); 00145 } 00146 } 00147 00148 /*****************************************************************/ 00149 /* */ 00150 /* probabilistic latent semantic analysis (pLSA) */ 00151 /* see T Hofmann, Probabilistic Latent Semantic */ 00152 /* Indexing for details */ 00153 /* */ 00154 /*****************************************************************/ 00155 00156 /** \brief Option object for the \ref pLSA algorithm. 00157 */ 00158 class PLSAOptions 00159 { 00160 public: 00161 /** Initialize all options with default values. 00162 */ 00163 PLSAOptions() 00164 : min_rel_gain(1e-4), 00165 max_iterations(50), 00166 normalized_component_weights(true) 00167 {} 00168 00169 /** Maximum number of iterations which is performed by the pLSA algorithm. 00170 00171 default: 50 00172 */ 00173 PLSAOptions & maximumNumberOfIterations(unsigned int n) 00174 { 00175 vigra_precondition(n >= 1, 00176 "PLSAOptions::maximumNumberOfIterations(): number must be a positive integer."); 00177 max_iterations = n; 00178 return *this; 00179 } 00180 00181 /** Minimum relative gain which is required for the algorithm to continue the iterations. 00182 00183 default: 1e-4 00184 */ 00185 PLSAOptions & minimumRelativeGain(double g) 00186 { 00187 vigra_precondition(g >= 0.0, 00188 "PLSAOptions::minimumRelativeGain(): number must be positive or zero."); 00189 min_rel_gain = g; 00190 return *this; 00191 } 00192 00193 /** Normalize the entries of the zv result array. 00194 00195 If true, the columns of zv sum to one. Otherwise, they have the same 00196 column sum as the original feature matrix. 00197 00198 default: true 00199 */ 00200 PLSAOptions & normalizedComponentWeights(bool v = true) 00201 { 00202 normalized_component_weights = v; 00203 return *this; 00204 } 00205 00206 double min_rel_gain; 00207 int max_iterations; 00208 bool normalized_component_weights; 00209 }; 00210 00211 /** \brief Decompose a matrix according to the pLSA algorithm. 00212 00213 This function implements the pLSA algorithm (probabilistic latent semantic analysis) 00214 proposed in 00215 00216 T. Hofmann: <a href="http://www.cs.brown.edu/people/th/papers/Hofmann-UAI99.pdf"> 00217 <i>"Probabilistic Latent Semantic Analysis"</i></a>, 00218 in: UAI'99, Proc. 15th Conf. on Uncertainty in Artificial Intelligence, 00219 pp. 289-296, Morgan Kaufmann, 1999 00220 00221 \arg features must be a matrix with shape <tt>(numFeatures * numSamples)</tt>, which is 00222 decomposed into the matrices 00223 \arg fz with shape <tt>(numFeatures * numComponents)</tt> and 00224 \arg zv with shape <tt>(numComponents * numSamples)</tt> 00225 00226 such that 00227 \f[ 00228 \mathrm{features} \approx \mathrm{fz} * \mathrm{zv} 00229 \f] 00230 (this formula applies when pLSA is called with 00231 <tt>PLSAOptions.normalizedComponentWeights(false)</tt>. Otherwise, you must 00232 normalize the features by calling <tt>\ref linalg::prepareColumns (features, features, UnitSum)</tt> 00233 to make the formula hold). 00234 00235 The shape parameter <tt>numComponents</tt> determines the complexity of 00236 the decomposition model and therefore the approximation quality. 00237 Intuitively, features are a set of words, and the samples a set of 00238 documents. The entries of the <tt>features</tt> matrix denote the relative 00239 frequency of the words in each document. The components represents a 00240 (presumably small) set of topics. The matrix <tt>fz</tt> encodes the 00241 relative frequency of words in the different topics, and the matrix 00242 <tt>zv</tt> encodes to what extend each topic explains the content of each 00243 document. 00244 00245 The option object determines the iteration termination conditions and the output 00246 normalization. In addition, you may pass a random number generator to pLSA() 00247 which is used to create the initial solution. 00248 00249 <b>Declarations:</b> 00250 00251 <b>\#include</b> <vigra/unsupervised_decomposition.hxx> 00252 00253 \code 00254 namespace vigra { 00255 00256 template <class U, class C1, class C2, class C3, class Random> 00257 void 00258 pLSA(MultiArrayView<2, U, C1> const & features, 00259 MultiArrayView<2, U, C2> & fz, 00260 MultiArrayView<2, U, C3> & zv, 00261 Random const& random, 00262 PLSAOptions const & options = PLSAOptions()); 00263 00264 template <class U, class C1, class C2, class C3> 00265 void 00266 pLSA(MultiArrayView<2, U, C1> const & features, 00267 MultiArrayView<2, U, C2> & fz, 00268 MultiArrayView<2, U, C3> & zv, 00269 PLSAOptions const & options = PLSAOptions()); 00270 } 00271 \endcode 00272 00273 <b>Usage:</b> 00274 \code 00275 Matrix<double> words(numWords, numDocuments); 00276 ... // fill the input matrix 00277 00278 int numTopics = 3; 00279 Matrix<double> fz(numWords, numTopics), 00280 zv(numTopics, numDocuments); 00281 00282 pLSA(words, fz, zv, PLSAOptions().normalizedComponentWeights(false)); 00283 00284 Matrix<double> model = fz*zv; 00285 double meanSquaredError = (words - model).squaredNorm() / numDocuments; 00286 \endcode 00287 */ 00288 doxygen_overloaded_function(template <...> void pLSA) 00289 00290 00291 template <class U, class C1, class C2, class C3, class Random> 00292 void 00293 pLSA(MultiArrayView<2, U, C1> const & features, 00294 MultiArrayView<2, U, C2> fz, 00295 MultiArrayView<2, U, C3> zv, 00296 Random const& random, 00297 PLSAOptions const & options = PLSAOptions()) 00298 { 00299 using namespace linalg; // activate matrix multiplication and arithmetic functions 00300 00301 int numFeatures = rowCount(features); 00302 int numSamples = columnCount(features); 00303 int numComponents = columnCount(fz); 00304 vigra_precondition(numFeatures >= numComponents && numComponents >= 1, 00305 "pLSA(): The number of features has to be larger or equal to the number of components in which the feature matrix is decomposed."); 00306 vigra_precondition(rowCount(fz) == numFeatures, 00307 "pLSA(): The output matrix fz has to be of dimension numFeatures*numComponents."); 00308 vigra_precondition(columnCount(zv) == numSamples && rowCount(zv) == numComponents, 00309 "pLSA(): The output matrix zv has to be of dimension numComponents*numSamples."); 00310 00311 // random initialization of result matrices, subsequent normalization 00312 UniformRandomFunctor<Random> randf(random); 00313 initMultiArray(destMultiArrayRange(fz), randf); 00314 initMultiArray(destMultiArrayRange(zv), randf); 00315 prepareColumns(fz, fz, UnitSum); 00316 prepareColumns(zv, zv, UnitSum); 00317 00318 // init vars 00319 double eps = 1.0/NumericTraits<U>::max(); // epsilon > 0 00320 double lastChange = NumericTraits<U>::max(); // infinity 00321 double err = 0; 00322 double err_old; 00323 int iteration = 0; 00324 00325 // expectation maximization (EM) algorithm 00326 Matrix<U> columnSums(1, numSamples); 00327 features.sum(columnSums); 00328 Matrix<U> expandedSums = ones<U>(numFeatures, 1) * columnSums; 00329 00330 while(iteration < options.max_iterations && (lastChange > options.min_rel_gain)) 00331 { 00332 Matrix<U> fzv = fz*zv; 00333 00334 //if(iteration%25 == 0) 00335 //{ 00336 //std::cout << "iteration: " << iteration << std::endl; 00337 //std::cout << "last relative change: " << lastChange << std::endl; 00338 //} 00339 00340 Matrix<U> factor = features / pointWise(fzv + (U)eps); 00341 zv *= (fz.transpose() * factor); 00342 fz *= (factor * zv.transpose()); 00343 prepareColumns(fz, fz, UnitSum); 00344 prepareColumns(zv, zv, UnitSum); 00345 00346 // check relative change in least squares model fit 00347 Matrix<U> model = expandedSums * pointWise(fzv); 00348 err_old = err; 00349 err = (features - model).squaredNorm(); 00350 //std::cout << "error: " << err << std::endl; 00351 lastChange = abs((err-err_old) / (U)(err + eps)); 00352 //std::cout << "lastChange: " << lastChange << std::endl; 00353 00354 iteration += 1; 00355 } 00356 //std::cout << "Terminated after " << iteration << " iterations." << std::endl; 00357 //std::cout << "Last relative change was " << lastChange << "." << std::endl; 00358 00359 if(!options.normalized_component_weights) 00360 { 00361 // undo the normalization 00362 for(int k=0; k<numSamples; ++k) 00363 columnVector(zv, k) *= columnSums(0, k); 00364 } 00365 } 00366 00367 template <class U, class C1, class C2, class C3> 00368 inline void 00369 pLSA(MultiArrayView<2, U, C1> const & features, 00370 MultiArrayView<2, U, C2> & fz, 00371 MultiArrayView<2, U, C3> & zv, 00372 PLSAOptions const & options = PLSAOptions()) 00373 { 00374 RandomNumberGenerator<> generator(RandomSeed); 00375 pLSA(features, fz, zv, generator, options); 00376 } 00377 00378 //@} 00379 00380 } // namespace vigra 00381 00382 00383 #endif // VIGRA_UNSUPERVISED_DECOMPOSITION_HXX 00384
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|