[ 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&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)