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

vigra/unsupervised_decomposition.hxx VIGRA

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&nbsp;(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> and 
00222         non-negative entries, which is 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&nbsp;(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)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.9.0 (Tue Nov 6 2012)