Tapkee
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
eigendecomposition.hpp
Go to the documentation of this file.
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
18 #endif
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 template <class MatrixType, class MatrixOperationType>
31 EigendecompositionResult eigendecomposition_impl_arpack(const MatrixType& wm, IndexType target_dimension, unsigned int skip)
32 {
33  timed_context context("ARPACK eigendecomposition");
34 
36  arpack(wm,target_dimension+skip,MatrixOperationType::ARPACK_CODE);
37 
38  if (arpack.info() == Eigen::Success)
39  {
40  std::stringstream ss;
41  ss << "Took " << arpack.getNbrIterations() << " iterations.";
43  DenseMatrix selected_eigenvectors = arpack.eigenvectors().rightCols(target_dimension);
44  return EigendecompositionResult(selected_eigenvectors,arpack.eigenvalues().tail(target_dimension));
45  }
46  else
47  {
48  throw eigendecomposition_error("eigendecomposition failed");
49  }
50  return EigendecompositionResult();
51 }
52 #endif
53 
55 template <class MatrixType, class MatrixOperationType>
56 EigendecompositionResult eigendecomposition_impl_dense(const MatrixType& wm, IndexType target_dimension, unsigned int skip)
57 {
58  timed_context context("Eigen library dense eigendecomposition");
59 
60  DenseMatrix dense_wm = wm;
61  DenseSelfAdjointEigenSolver solver(dense_wm);
62 
63  if (solver.info() == Eigen::Success)
64  {
65  if (MatrixOperationType::largest)
66  {
67  assert(skip==0);
68  DenseMatrix selected_eigenvectors = solver.eigenvectors().rightCols(target_dimension);
69  return EigendecompositionResult(selected_eigenvectors,solver.eigenvalues().tail(target_dimension));
70  }
71  else
72  {
73  DenseMatrix selected_eigenvectors = solver.eigenvectors().leftCols(target_dimension+skip).rightCols(target_dimension);
74  return EigendecompositionResult(selected_eigenvectors,solver.eigenvalues().segment(skip,skip+target_dimension));
75  }
76  }
77  else
78  {
79  throw eigendecomposition_error("eigendecomposition failed");
80  }
81  return EigendecompositionResult();
82 }
83 
85 template <class MatrixType, class MatrixOperationType>
86 EigendecompositionResult eigendecomposition_impl_randomized(const MatrixType& wm, IndexType target_dimension, unsigned int skip)
87 {
88  timed_context context("Randomized eigendecomposition");
89 
90  DenseMatrix O(wm.rows(), target_dimension+skip);
91  for (IndexType i=0; i<O.rows(); ++i)
92  {
93  for (IndexType j=0; j<O.cols(); j++)
94  {
95  O(i,j) = tapkee::gaussian_random();
96  }
97  }
98  MatrixOperationType operation(wm);
99 
100  DenseMatrix Y = operation(O);
101  for (IndexType i=0; i<Y.cols(); i++)
102  {
103  for (IndexType j=0; j<i; j++)
104  {
105  ScalarType r = Y.col(i).dot(Y.col(j));
106  Y.col(i) -= r*Y.col(j);
107  }
108  ScalarType norm = Y.col(i).norm();
109  if (norm < 1e-4)
110  {
111  for (int k = i; k<Y.cols(); k++)
112  Y.col(k).setZero();
113  }
114  Y.col(i) *= (1.f / norm);
115  }
116 
117  DenseMatrix B1 = operation(Y);
118  DenseMatrix B = Y.householderQr().solve(B1);
119  DenseSelfAdjointEigenSolver eigenOfB(B);
120 
121  if (eigenOfB.info() == Eigen::Success)
122  {
123  if (MatrixOperationType::largest)
124  {
125  assert(skip==0);
126  DenseMatrix selected_eigenvectors = (Y*eigenOfB.eigenvectors()).rightCols(target_dimension);
127  return EigendecompositionResult(selected_eigenvectors,eigenOfB.eigenvalues());
128  }
129  else
130  {
131  DenseMatrix selected_eigenvectors = (Y*eigenOfB.eigenvectors()).leftCols(target_dimension+skip).rightCols(target_dimension);
132  return EigendecompositionResult(selected_eigenvectors,eigenOfB.eigenvalues());
133  }
134  }
135  else
136  {
137  throw eigendecomposition_error("eigendecomposition failed");
138  }
139  return EigendecompositionResult();
140 }
141 
170 template <class MatrixType, class MatrixOperationType>
172  IndexType target_dimension, unsigned int skip)
173 {
174  LoggingSingleton::instance().message_info("Using the " + get_eigen_method_name(method) + " eigendecomposition method.");
175  switch (method)
176  {
177 #ifdef TAPKEE_WITH_ARPACK
178  case Arpack:
179  return eigendecomposition_impl_arpack<MatrixType, MatrixOperationType>(m, target_dimension, skip);
180 #endif
181  case Randomized:
182  return eigendecomposition_impl_randomized<MatrixType, MatrixOperationType>(m, target_dimension, skip);
183  case Dense:
184  return eigendecomposition_impl_dense<MatrixType, MatrixOperationType>(m, target_dimension, skip);
185  default: break;
186  }
187  return EigendecompositionResult();
188 }
189 
190 } // End of namespace tapkee_internal
191 } // End of namespace tapkee
192 
193 #endif
EigendecompositionResult eigendecomposition_impl_randomized(const MatrixType &wm, IndexType target_dimension, unsigned int skip)
Randomized redsvd-like implementation of eigendecomposition-based embedding.
EigendecompositionResult eigendecomposition(EigenMethod method, const MatrixType &m, IndexType target_dimension, unsigned int skip)
Multiple implementation handler method for various eigendecomposition methods.
Eigen::Matrix< tapkee::ScalarType, Eigen::Dynamic, Eigen::Dynamic > DenseMatrix
dense matrix type (non-overridable)
Definition: types.hpp:23
ARPACK-based method (requires the ARPACK library binaries to be available around). Recommended to be used as a default method. Supports both generalized and standard eigenproblems.
Eigen library dense method (could be useful for debugging). Computes all eigenvectors thus can be ver...
std::string get_eigen_method_name(EigenMethod m)
Definition: naming.hpp:56
double ScalarType
default scalar value (can be overrided with TAPKEE_CUSTOM_INTERNAL_NUMTYPE define) ...
Definition: types.hpp:15
An exception type that is thrown when eigendecomposition is failed.
Definition: exceptions.hpp:85
Randomized method (implementation taken from the redsvd lib). Supports only standard but not generali...
EigenMethod
Eigendecomposition methods.
int IndexType
indexing type (non-overridable) set to int for compatibility with OpenMP 2.0
Definition: types.hpp:19
EigendecompositionResult eigendecomposition_impl_dense(const MatrixType &wm, IndexType target_dimension, unsigned int skip)
Eigen library dense implementation of eigendecomposition-based embedding.
Eigen::SelfAdjointEigenSolver< tapkee::DenseMatrix > DenseSelfAdjointEigenSolver
selfadjoint solver (non-overridable)
Definition: types.hpp:33
const Matrix< Scalar, Dynamic, Dynamic > & eigenvectors() const
Returns the eigenvectors of given matrix.
void message_info(const std::string &msg)
Definition: logging.hpp:113
const Matrix< Scalar, Dynamic, 1 > & eigenvalues() const
Returns the eigenvalues of given matrix.
ComputationInfo info() const
Reports whether previous computation was successful.
static LoggingSingleton & instance()
Definition: logging.hpp:100
EigendecompositionResult eigendecomposition_impl_arpack(const MatrixType &wm, IndexType target_dimension, unsigned int skip)
ARPACK implementation of eigendecomposition-based embedding.
ScalarType gaussian_random()
Definition: random.hpp:39
TAPKEE_INTERNAL_PAIR< tapkee::DenseMatrix, tapkee::DenseVector > EigendecompositionResult
Definition: synonyms.hpp:41