LCOV - code coverage report
Current view: top level - routines - eigendecomposition.hpp (source / functions) Hit Total Coverage
Test: clean.info Lines: 54 59 91.5 %
Date: 2013-05-24 Functions: 11 16 68.8 %
Branches: 224 864 25.9 %

           Branch data     Line data    Source code
       1                 :            : /* This software is distributed under BSD 3-clause license (see LICENSE file).
       2                 :            :  *
       3                 :            :  * Copyright (c) 2012-2013 Sergey Lisitsyn
       4                 :            :  *
       5                 :            :  * Randomized eigendecomposition code is inspired by the redsvd library
       6                 :            :  * code which is also distributed under BSD 3-clause license.
       7                 :            :  *
       8                 :            :  * Copyright (c) 2010-2013 Daisuke Okanohara
       9                 :            :  *
      10                 :            :  */
      11                 :            : 
      12                 :            : #ifndef TAPKEE_EIGENDECOMPOSITION_H_
      13                 :            : #define TAPKEE_EIGENDECOMPOSITION_H_
      14                 :            : 
      15                 :            : /* Tapkee includes */
      16                 :            : #ifdef TAPKEE_WITH_ARPACK
      17                 :            :         #include <tapkee/utils/arpack_wrapper.hpp>
      18                 :            : #endif
      19                 :            : #include <tapkee/routines/matrix_operations.hpp>
      20                 :            : #include <tapkee/defines.hpp>
      21                 :            : /* End of Tapkee includes */
      22                 :            : 
      23                 :            : namespace tapkee
      24                 :            : {
      25                 :            : namespace tapkee_internal
      26                 :            : {
      27                 :            : 
      28                 :            : #ifdef TAPKEE_WITH_ARPACK
      29                 :            : //! ARPACK implementation of eigendecomposition-based embedding
      30                 :            : template <class MatrixType, class MatrixOperationType> 
      31                 :         23 : EigendecompositionResult eigendecomposition_impl_arpack(const MatrixType& wm, IndexType target_dimension, unsigned int skip)
      32                 :            : {
      33 [ +  - ][ +  - ]:         23 :         timed_context context("ARPACK eigendecomposition");
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
      34                 :            : 
      35                 :            :         ArpackGeneralizedSelfAdjointEigenSolver<MatrixType, MatrixType, MatrixOperationType> 
      36 [ +  - ][ +  - ]:         23 :                 arpack(wm,target_dimension+skip,MatrixOperationType::ARPACK_CODE);
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
      37                 :            : 
      38   [ +  -  +  -  :         23 :         if (arpack.info() == Eigen::Success)
             +  -  +  - ]
      39                 :            :         {
      40 [ +  - ][ +  - ]:         23 :                 std::stringstream ss;
         [ +  - ][ +  - ]
      41 [ +  - ][ +  - ]:         23 :                 ss << "Took " << arpack.getNbrIterations() << " iterations.";
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
      42 [ +  - ][ +  - ]:         23 :                 LoggingSingleton::instance().message_info(ss.str());
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
      43 [ +  - ][ +  - ]:         23 :                 DenseMatrix selected_eigenvectors = arpack.eigenvectors().rightCols(target_dimension);
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
      44 [ +  - ][ +  - ]:         23 :                 return EigendecompositionResult(selected_eigenvectors,arpack.eigenvalues().tail(target_dimension));
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
      45                 :            :         }
      46                 :            :         else
      47                 :            :         {
      48 [ #  # ][ #  # ]:          0 :                 throw eigendecomposition_error("eigendecomposition failed");
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
      49                 :            :         }
      50 [ +  - ][ +  - ]:         23 :         return EigendecompositionResult();
         [ +  - ][ +  - ]
      51                 :            : }
      52                 :            : #endif
      53                 :            : 
      54                 :            : //! Eigen library dense implementation of eigendecomposition-based embedding
      55                 :            : template <class MatrixType, class MatrixOperationType> 
      56                 :          2 : EigendecompositionResult eigendecomposition_impl_dense(const MatrixType& wm, IndexType target_dimension, unsigned int skip)
      57                 :            : {
      58 [ +  - ][ +  - ]:          2 :         timed_context context("Eigen library dense eigendecomposition");
         [ +  - ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ +  - ][ +  - ]
                 [ +  - ]
      59                 :            : 
      60   [ +  -  +  -  :          2 :         DenseMatrix dense_wm = wm;
             #  #  #  # ]
      61 [ +  - ][ +  - ]:          2 :         DenseSelfAdjointEigenSolver solver(dense_wm);
         [ #  # ][ #  # ]
      62                 :            : 
      63 [ +  - ][ +  - ]:          2 :         if (solver.info() == Eigen::Success)
         [ #  # ][ #  # ]
      64                 :            :         {
      65                 :            :                 if (MatrixOperationType::largest)
      66                 :            :                 {
      67 [ -  + ][ #  # ]:          1 :                         assert(skip==0);
                 [ #  # ]
      68 [ +  - ][ +  - ]:          1 :                         DenseMatrix selected_eigenvectors = solver.eigenvectors().rightCols(target_dimension);
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
      69 [ +  - ][ +  - ]:          1 :                         return EigendecompositionResult(selected_eigenvectors,solver.eigenvalues().tail(target_dimension));
         [ +  - ][ +  - ]
         [ +  - ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ #  # ]
      70                 :            :                 } 
      71                 :            :                 else
      72                 :            :                 {
      73 [ +  - ][ +  - ]:          1 :                         DenseMatrix selected_eigenvectors = solver.eigenvectors().leftCols(target_dimension+skip).rightCols(target_dimension);
                 [ +  - ]
      74 [ +  - ][ +  - ]:          1 :                         return EigendecompositionResult(selected_eigenvectors,solver.eigenvalues().segment(skip,skip+target_dimension));
         [ +  - ][ +  - ]
                 [ +  - ]
      75                 :            :                 }
      76                 :            :         }
      77                 :            :         else
      78                 :            :         {
      79 [ #  # ][ #  # ]:          0 :                 throw eigendecomposition_error("eigendecomposition failed");
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
      80                 :            :         }
      81 [ +  - ][ +  - ]:          2 :         return EigendecompositionResult();
         [ +  - ][ +  - ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
      82                 :            : }
      83                 :            : 
      84                 :            : //! Randomized redsvd-like implementation of eigendecomposition-based embedding
      85                 :            : template <class MatrixType, class MatrixOperationType> 
      86                 :          2 : EigendecompositionResult eigendecomposition_impl_randomized(const MatrixType& wm, IndexType target_dimension, unsigned int skip)
      87                 :            : {
      88 [ #  # ][ #  # ]:          2 :         timed_context context("Randomized eigendecomposition");
         [ #  # ][ +  - ]
         [ +  - ][ +  - ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
      89                 :            :         
      90         [ #  # ]:          2 :         DenseMatrix O(wm.rows(), target_dimension+skip);
           [ #  #  +  + ]
           [ +  -  #  # ]
           [ #  #  #  # ]
         [ #  +  - ][ - ]
      91 [ #  # ][ #  # ]:        104 :         for (IndexType i=0; i<O.rows(); ++i)
         [ +  - ][ +  + ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
      92                 :            :         {
      93 [ #  # ][ #  # ]:        304 :                 for (IndexType j=0; j<O.cols(); j++)
         [ +  - ][ +  + ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
      94                 :            :                 {
      95 [ #  # ][ +  - ]:        202 :                         O(i,j) = tapkee::gaussian_random();
         [ #  # ][ #  # ]
      96                 :            :                 }
      97                 :            :         }
      98 [ #  # ][ #  # ]:          2 :         MatrixOperationType operation(wm);
      99                 :            : 
     100   [ #  #  +  -  :          2 :         DenseMatrix Y = operation(O);
           #  # ][ #  #  
             #  #  +  - ]
     101 [ #  # ][ #  # ]:          5 :         for (IndexType i=0; i<Y.cols(); i++)
         [ +  - ][ +  + ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
     102                 :            :         {
     103 [ #  # ][ +  + ]:          4 :                 for (IndexType j=0; j<i; j++)
         [ #  # ][ #  # ]
     104                 :            :                 {
     105 [ #  # ][ #  # ]:          1 :                         ScalarType r = Y.col(i).dot(Y.col(j));
         [ #  # ][ +  - ]
         [ +  - ][ +  - ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
     106 [ #  # ][ #  # ]:          1 :                         Y.col(i) -= r*Y.col(j);
         [ #  # ][ #  # ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
     107                 :            :                 }
     108 [ #  # ][ #  # ]:          3 :                 ScalarType norm = Y.col(i).norm();
         [ +  - ][ +  - ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
     109 [ #  # ][ +  + ]:          3 :                 if (norm < 1e-4)
         [ #  # ][ #  # ]
     110                 :            :                 {
     111 [ #  # ][ #  # ]:          3 :                         for (int k = i; k<Y.cols(); k++)
         [ +  - ][ +  + ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
     112 [ #  # ][ #  # ]:          2 :                                 Y.col(k).setZero();
         [ +  - ][ +  - ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
     113                 :            :                 }
     114 [ #  # ][ #  # ]:          3 :                 Y.col(i) *= (1.f / norm);
         [ +  - ][ +  - ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
     115                 :            :         }
     116                 :            : 
     117 [ #  # ][ +  - ]:          2 :         DenseMatrix B1 = operation(Y);
         [ #  # ][ #  # ]
     118 [ #  # ][ #  # ]:          2 :         DenseMatrix B = Y.householderQr().solve(B1);
         [ #  # ][ #  # ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
     119 [ #  # ][ +  - ]:          2 :         DenseSelfAdjointEigenSolver eigenOfB(B);
         [ #  # ][ #  # ]
     120                 :            : 
     121 [ #  # ][ +  + ]:          2 :         if (eigenOfB.info() == Eigen::Success)
         [ #  # ][ #  # ]
     122                 :            :         {
     123                 :            :                 if (MatrixOperationType::largest)
     124                 :            :                 {
     125 [ -  + ][ #  # ]:          1 :                         assert(skip==0);
                 [ #  # ]
     126 [ +  - ][ +  - ]:          1 :                         DenseMatrix selected_eigenvectors = (Y*eigenOfB.eigenvectors()).rightCols(target_dimension);
         [ +  - ][ +  - ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
     127 [ +  - ][ +  - ]:          1 :                         return EigendecompositionResult(selected_eigenvectors,eigenOfB.eigenvalues());
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
     128                 :            :                 } 
     129                 :            :                 else
     130                 :            :                 {
     131 [ #  # ][ #  # ]:          0 :                         DenseMatrix selected_eigenvectors = (Y*eigenOfB.eigenvectors()).leftCols(target_dimension+skip).rightCols(target_dimension);
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
     132 [ #  # ][ #  # ]:          0 :                         return EigendecompositionResult(selected_eigenvectors,eigenOfB.eigenvalues());
     133                 :            :                 }
     134                 :            :         }
     135                 :            :         else
     136                 :            :         {
     137 [ #  # ][ #  # ]:          1 :                 throw eigendecomposition_error("eigendecomposition failed");
         [ +  - ][ +  - ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
     138                 :            :         }
     139 [ #  # ][ #  # ]:          2 :         return EigendecompositionResult();
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
                 [ +  - ]
     140                 :            : }
     141                 :            : 
     142                 :            : //! Multiple implementation handler method for various eigendecomposition methods. 
     143                 :            : //!
     144                 :            : //! Has three template parameters:
     145                 :            : //! MatrixType - class of weight matrix to perform eigendecomposition of
     146                 :            : //! MatrixOperationType - class of product operation over matrix.
     147                 :            : //!
     148                 :            : //! In order to compute largest eigenvalues MatrixOperationType should provide
     149                 :            : //! implementation of operator()(DenseMatrix) which computes right product
     150                 :            : //! of the parameter with the MatrixType.
     151                 :            : //! 
     152                 :            : //! In order to compute smallest eigenvalues MatrixOperationType should provide
     153                 :            : //! implementation of operator()(DenseMatrix) which solves linear system with
     154                 :            : //! given right-hand side part. 
     155                 :            : //! 
     156                 :            : //! Currently supports three methods:
     157                 :            : //!
     158                 :            : //! <ul>
     159                 :            : //! <li> Arpack
     160                 :            : //! <li> Randomized
     161                 :            : //! <li> Dense
     162                 :            : //! </ul>
     163                 :            : //!
     164                 :            : //! @param method one of supported eigendecomposition methods
     165                 :            : //! @param m matrix to be eigendecomposed 
     166                 :            : //! @param target_dimension target dimension of embedding i.e. number of eigenvectors to be
     167                 :            : //!        computed
     168                 :            : //! @param skip number of eigenvectors to skip (from either smallest or largest side)
     169                 :            : //!
     170                 :            : template <class MatrixType, class MatrixOperationType>
     171                 :         27 : EigendecompositionResult eigendecomposition(EigenMethod method, const MatrixType& m, 
     172                 :            :                                             IndexType target_dimension, unsigned int skip)
     173                 :            : {
     174 [ +  - ][ +  - ]:         27 :         LoggingSingleton::instance().message_info("Using the " + get_eigen_method_name(method) + " eigendecomposition method.");
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
     175   [ +  -  +  -  :         27 :         switch (method)
          +  +  +  -  +  
          -  -  -  +  -  
                   -  - ]
     176                 :            :         {
     177                 :            : #ifdef TAPKEE_WITH_ARPACK
     178                 :            :                 case Arpack: 
     179                 :         23 :                         return eigendecomposition_impl_arpack<MatrixType, MatrixOperationType>(m, target_dimension, skip);
     180                 :            : #endif
     181                 :            :                 case Randomized: 
     182                 :          2 :                         return eigendecomposition_impl_randomized<MatrixType, MatrixOperationType>(m, target_dimension, skip);
     183                 :            :                 case Dense:
     184                 :          2 :                         return eigendecomposition_impl_dense<MatrixType, MatrixOperationType>(m, target_dimension, skip);
     185                 :          0 :                 default: break;
     186                 :            :         }
     187                 :         26 :         return EigendecompositionResult();
     188                 :            : }
     189                 :            : 
     190                 :            : } // End of namespace tapkee_internal
     191                 :            : } // End of namespace tapkee
     192                 :            : 
     193                 :            : #endif

Generated by: LCOV version 1.9