Tapkee
eigendecomposition.hpp
Go to the documentation of this file.
00001 /* This software is distributed under BSD 3-clause license (see LICENSE file).
00002  *
00003  * Copyright (c) 2012-2013 Sergey Lisitsyn
00004  *
00005  * Randomized eigendecomposition code is inspired by the redsvd library
00006  * code which is also distributed under BSD 3-clause license.
00007  *
00008  * Copyright (c) 2010-2013 Daisuke Okanohara
00009  *
00010  */
00011 
00012 #ifndef TAPKEE_EIGENDECOMPOSITION_H_
00013 #define TAPKEE_EIGENDECOMPOSITION_H_
00014 
00015 /* Tapkee includes */
00016 #ifdef TAPKEE_WITH_ARPACK
00017     #include <tapkee/utils/arpack_wrapper.hpp>
00018 #endif
00019 #include <tapkee/routines/matrix_operations.hpp>
00020 #include <tapkee/defines.hpp>
00021 /* End of Tapkee includes */
00022 
00023 namespace tapkee
00024 {
00025 namespace tapkee_internal
00026 {
00027 
00028 #ifdef TAPKEE_WITH_ARPACK
00029 
00030 template <class MatrixType, class MatrixOperationType> 
00031 EigendecompositionResult eigendecomposition_impl_arpack(const MatrixType& wm, IndexType target_dimension, unsigned int skip)
00032 {
00033     timed_context context("ARPACK eigendecomposition");
00034 
00035     ArpackGeneralizedSelfAdjointEigenSolver<MatrixType, MatrixType, MatrixOperationType> 
00036         arpack(wm,target_dimension+skip,MatrixOperationType::ARPACK_CODE);
00037 
00038     if (arpack.info() == Eigen::Success)
00039     {
00040         std::stringstream ss;
00041         ss << "Took " << arpack.getNbrIterations() << " iterations.";
00042         LoggingSingleton::instance().message_info(ss.str());
00043         DenseMatrix selected_eigenvectors = arpack.eigenvectors().rightCols(target_dimension);
00044         return EigendecompositionResult(selected_eigenvectors,arpack.eigenvalues().tail(target_dimension));
00045     }
00046     else
00047     {
00048         throw eigendecomposition_error("eigendecomposition failed");
00049     }
00050     return EigendecompositionResult();
00051 }
00052 #endif
00053 
00055 template <class MatrixType, class MatrixOperationType> 
00056 EigendecompositionResult eigendecomposition_impl_dense(const MatrixType& wm, IndexType target_dimension, unsigned int skip)
00057 {
00058     timed_context context("Eigen library dense eigendecomposition");
00059 
00060     DenseMatrix dense_wm = wm;
00061     DenseSelfAdjointEigenSolver solver(dense_wm);
00062 
00063     if (solver.info() == Eigen::Success)
00064     {
00065         if (MatrixOperationType::largest)
00066         {
00067             assert(skip==0);
00068             DenseMatrix selected_eigenvectors = solver.eigenvectors().rightCols(target_dimension);
00069             return EigendecompositionResult(selected_eigenvectors,solver.eigenvalues().tail(target_dimension));
00070         } 
00071         else
00072         {
00073             DenseMatrix selected_eigenvectors = solver.eigenvectors().leftCols(target_dimension+skip).rightCols(target_dimension);
00074             return EigendecompositionResult(selected_eigenvectors,solver.eigenvalues().segment(skip,skip+target_dimension));
00075         }
00076     }
00077     else
00078     {
00079         throw eigendecomposition_error("eigendecomposition failed");
00080     }
00081     return EigendecompositionResult();
00082 }
00083 
00085 template <class MatrixType, class MatrixOperationType> 
00086 EigendecompositionResult eigendecomposition_impl_randomized(const MatrixType& wm, IndexType target_dimension, unsigned int skip)
00087 {
00088     timed_context context("Randomized eigendecomposition");
00089     
00090     DenseMatrix O(wm.rows(), target_dimension+skip);
00091     for (IndexType i=0; i<O.rows(); ++i)
00092     {
00093         for (IndexType j=0; j<O.cols(); j++)
00094         {
00095             O(i,j) = tapkee::gaussian_random();
00096         }
00097     }
00098     MatrixOperationType operation(wm);
00099 
00100     DenseMatrix Y = operation(O);
00101     for (IndexType i=0; i<Y.cols(); i++)
00102     {
00103         for (IndexType j=0; j<i; j++)
00104         {
00105             ScalarType r = Y.col(i).dot(Y.col(j));
00106             Y.col(i) -= r*Y.col(j);
00107         }
00108         ScalarType norm = Y.col(i).norm();
00109         if (norm < 1e-4)
00110         {
00111             for (int k = i; k<Y.cols(); k++)
00112                 Y.col(k).setZero();
00113         }
00114         Y.col(i) *= (1.f / norm);
00115     }
00116 
00117     DenseMatrix B1 = operation(Y);
00118     DenseMatrix B = Y.householderQr().solve(B1);
00119     DenseSelfAdjointEigenSolver eigenOfB(B);
00120 
00121     if (eigenOfB.info() == Eigen::Success)
00122     {
00123         if (MatrixOperationType::largest)
00124         {
00125             assert(skip==0);
00126             DenseMatrix selected_eigenvectors = (Y*eigenOfB.eigenvectors()).rightCols(target_dimension);
00127             return EigendecompositionResult(selected_eigenvectors,eigenOfB.eigenvalues());
00128         } 
00129         else
00130         {
00131             DenseMatrix selected_eigenvectors = (Y*eigenOfB.eigenvectors()).leftCols(target_dimension+skip).rightCols(target_dimension);
00132             return EigendecompositionResult(selected_eigenvectors,eigenOfB.eigenvalues());
00133         }
00134     }
00135     else
00136     {
00137         throw eigendecomposition_error("eigendecomposition failed");
00138     }
00139     return EigendecompositionResult();
00140 }
00141 
00170 template <class MatrixType, class MatrixOperationType>
00171 EigendecompositionResult eigendecomposition(EigenMethod method, const MatrixType& m, 
00172                                             IndexType target_dimension, unsigned int skip)
00173 {
00174     LoggingSingleton::instance().message_info("Using the " + get_eigen_method_name(method) + " eigendecomposition method.");
00175     switch (method)
00176     {
00177 #ifdef TAPKEE_WITH_ARPACK
00178         case Arpack: 
00179             return eigendecomposition_impl_arpack<MatrixType, MatrixOperationType>(m, target_dimension, skip);
00180 #endif
00181         case Randomized: 
00182             return eigendecomposition_impl_randomized<MatrixType, MatrixOperationType>(m, target_dimension, skip);
00183         case Dense:
00184             return eigendecomposition_impl_dense<MatrixType, MatrixOperationType>(m, target_dimension, skip);
00185         default: break;
00186     }
00187     return EigendecompositionResult();
00188 }
00189 
00190 } // End of namespace tapkee_internal
00191 } // End of namespace tapkee
00192 
00193 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines