Tapkee
generalized_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 
6 #ifndef TAPKEE_GENERALIZED_EIGENDECOMPOSITION_H_
7 #define TAPKEE_GENERALIZED_EIGENDECOMPOSITION_H_
8 
9 /* Tapkee includes */
10 #ifdef TAPKEE_WITH_ARPACK
12 #endif
14 /* End of Tapkee includes */
15 
16 namespace tapkee
17 {
18 namespace tapkee_internal
19 {
20 
21 #ifdef TAPKEE_WITH_ARPACK
22 template <class LMatrixType, class RMatrixType, class MatrixOperationType>
25  const RMatrixType& rhs, IndexType target_dimension, unsigned int skip)
26 {
27  timed_context context("ARPACK DSXUPD generalized eigendecomposition");
28 
30  arpack(lhs,rhs,target_dimension+skip,"SM");
31 
32  if (arpack.info() == Eigen::Success)
33  {
34  std::string message = formatting::format("Took {} iterations.", arpack.getNbrIterations());
36  DenseMatrix selected_eigenvectors = (arpack.eigenvectors()).rightCols(target_dimension);
37  return EigendecompositionResult(selected_eigenvectors,arpack.eigenvalues().tail(target_dimension));
38  }
39  else
40  {
41  throw eigendecomposition_error("eigendecomposition failed");
42  }
43  return EigendecompositionResult();
44 }
45 #endif
46 
48 template <class LMatrixType, class RMatrixType, class MatrixOperationType>
50  const RMatrixType& rhs, IndexType target_dimension, unsigned int skip)
51 {
52  timed_context context("Eigen dense generalized eigendecomposition");
53 
54  DenseMatrix dense_lhs = lhs;
55  DenseMatrix dense_rhs = rhs;
56  Eigen::GeneralizedSelfAdjointEigenSolver<DenseMatrix> solver(dense_lhs, dense_rhs);
57  if (solver.info() == Eigen::Success)
58  {
59  if (MatrixOperationType::largest)
60  {
61  assert(skip==0);
62  DenseMatrix selected_eigenvectors = solver.eigenvectors().rightCols(target_dimension);
63  return EigendecompositionResult(selected_eigenvectors,solver.eigenvalues().tail(target_dimension));
64  }
65  else
66  {
67  DenseMatrix selected_eigenvectors = solver.eigenvectors().leftCols(target_dimension+skip).rightCols(target_dimension);
68  return EigendecompositionResult(selected_eigenvectors,solver.eigenvalues().segment(skip,skip+target_dimension));
69  }
70  }
71  else
72  {
73  throw eigendecomposition_error("eigendecomposition failed");
74  }
75 
76  return EigendecompositionResult();
77 }
78 
79 template <typename LMatrixType, typename RMatrixType>
81 {
82 #ifdef TAPKEE_WITH_ARPACK
83  EigendecompositionResult arpack(const LMatrixType& lhs, const RMatrixType& rhs,
84  const ComputationStrategy& strategy,
85  const EigendecompositionStrategy& eigen_strategy,
86  IndexType target_dimension);
87 #endif
88  EigendecompositionResult dense(const LMatrixType& lhs, const RMatrixType& rhs,
89  const ComputationStrategy& strategy,
90  const EigendecompositionStrategy& eigen_strategy,
91  IndexType target_dimension);
92 };
93 
94 template <>
96 {
97 #ifdef TAPKEE_WITH_ARPACK
99  const ComputationStrategy& strategy,
100  const EigendecompositionStrategy& eigen_strategy,
101  IndexType target_dimension)
102  {
103  if (strategy.is(HomogeneousCPUStrategy))
104  {
105  if (eigen_strategy.is(SmallestEigenvalues))
108  (lhs,rhs,target_dimension,eigen_strategy.skip());
109  unsupported();
110  }
111  unsupported();
112  return EigendecompositionResult();
113  }
114 #endif
116  const ComputationStrategy& strategy,
117  const EigendecompositionStrategy& eigen_strategy,
118  IndexType target_dimension)
119  {
120  if (strategy.is(HomogeneousCPUStrategy))
121  {
122  if (eigen_strategy.is(SmallestEigenvalues))
125  (lhs,rhs,target_dimension,eigen_strategy.skip());
126  unsupported();
127  }
128  unsupported();
129  return EigendecompositionResult();
130  }
131  inline void unsupported() const
132  {
133  throw unsupported_method_error("Unsupported method");
134  }
135 };
136 
137 template <>
139 {
140 #ifdef TAPKEE_WITH_ARPACK
142  const ComputationStrategy& strategy,
143  const EigendecompositionStrategy& eigen_strategy,
144  IndexType target_dimension)
145  {
146  if (strategy.is(HomogeneousCPUStrategy))
147  {
148  if (eigen_strategy.is(SmallestEigenvalues))
151  (lhs,rhs,target_dimension,0);
152  unsupported();
153  }
154  unsupported();
155  return EigendecompositionResult();
156  }
157 #endif
159  const ComputationStrategy& strategy,
160  const EigendecompositionStrategy& eigen_strategy,
161  IndexType target_dimension)
162  {
163  if (strategy.is(HomogeneousCPUStrategy))
164  {
165  if (eigen_strategy.is(SmallestEigenvalues))
168  (lhs,rhs,target_dimension,0);
169  unsupported();
170  }
171  unsupported();
172  return EigendecompositionResult();
173  }
174  inline void unsupported() const
175  {
176  throw unsupported_method_error("Unsupported method");
177  }
178 };
179 
180 template <class LMatrixType, class RMatrixType>
182  const EigendecompositionStrategy& eigen_strategy,
183  const LMatrixType& lhs, const RMatrixType& rhs, IndexType target_dimension)
184 {
185  LoggingSingleton::instance().message_info(formatting::format("Using the {} eigendecomposition method.",
186  get_eigen_method_name(method)));
187 #ifdef TAPKEE_WITH_ARPACK
188  if (method.is(Arpack))
190  .arpack(lhs, rhs, strategy, eigen_strategy, target_dimension);
191 #endif
192  if (method.is(Dense))
194  .dense(lhs, rhs, strategy, eigen_strategy, target_dimension);
195  if (method.is(Randomized))
196  throw unsupported_method_error("Randomized method is not supported for generalized eigenproblems");
197  return EigendecompositionResult();
198 }
199 
200 } // End of namespace tapkee_internal
201 } // End of namespace tapkee
202 
203 #endif
Eigen::DiagonalMatrix< tapkee::ScalarType, Eigen::Dynamic > DenseDiagonalMatrix
dense diagonal matrix
Definition: types.hpp:27
static const EigenMethod Dense("Dense")
Eigen library dense method (could be useful for debugging). Computes all eigenvectors thus can be ver...
static const EigendecompositionStrategy SmallestEigenvalues("Smallest eigenvalues", 1)
static const EigenMethod Randomized("Randomized")
Randomized method (implementation taken from the redsvd lib). Supports only standard but not generali...
Matrix-matrix operation used to compute smallest eigenvalues and associated eigenvectors of a dense m...
const Matrix< Scalar, Dynamic, Dynamic > & eigenvectors() const
Returns the eigenvectors of given matrix.
EigendecompositionResult dense(const SparseWeightMatrix &lhs, const DenseDiagonalMatrix &rhs, const ComputationStrategy &strategy, const EigendecompositionStrategy &eigen_strategy, IndexType target_dimension)
Eigen::Matrix< tapkee::ScalarType, Eigen::Dynamic, Eigen::Dynamic > DenseMatrix
dense matrix type (non-overridable)
Definition: types.hpp:23
EigendecompositionResult arpack(const LMatrixType &lhs, const RMatrixType &rhs, const ComputationStrategy &strategy, const EigendecompositionStrategy &eigen_strategy, IndexType target_dimension)
EigendecompositionResult dense(const LMatrixType &lhs, const RMatrixType &rhs, const ComputationStrategy &strategy, const EigendecompositionStrategy &eigen_strategy, IndexType target_dimension)
EigendecompositionResult generalized_eigendecomposition(const EigenMethod &method, const ComputationStrategy &strategy, const EigendecompositionStrategy &eigen_strategy, const LMatrixType &lhs, const RMatrixType &rhs, IndexType target_dimension)
ComputationInfo info() const
Reports whether previous computation was successful.
EigendecompositionResult arpack(const SparseWeightMatrix &lhs, const DenseDiagonalMatrix &rhs, const ComputationStrategy &strategy, const EigendecompositionStrategy &eigen_strategy, IndexType target_dimension)
An exception type that is thrown when eigendecomposition is failed.
EigendecompositionResult dense(const DenseMatrix &lhs, const DenseMatrix &rhs, const ComputationStrategy &strategy, const EigendecompositionStrategy &eigen_strategy, IndexType target_dimension)
Eigen::SparseMatrix< tapkee::ScalarType > SparseWeightMatrix
sparse weight matrix type (non-overridable)
Definition: types.hpp:29
static const EigenMethod Arpack("Arpack")
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.
std::string format(const std::string &fmt, const ValueWrapper a)
Definition: formatting.hpp:411
EigendecompositionResult generalized_eigendecomposition_impl_dense(const LMatrixType &lhs, const RMatrixType &rhs, IndexType target_dimension, unsigned int skip)
Eigen library dense implementation of eigendecomposition.
int IndexType
indexing type (non-overridable) set to int for compatibility with OpenMP 2.0
Definition: types.hpp:19
Matrix-matrix operation used to compute smallest eigenvalues and associated eigenvectors of a sparse ...
EigendecompositionResult arpack(const DenseMatrix &lhs, const DenseMatrix &rhs, const ComputationStrategy &strategy, const EigendecompositionStrategy &eigen_strategy, IndexType target_dimension)
static const ComputationStrategy HomogeneousCPUStrategy("CPU")
void message_info(const std::string &msg)
Definition: logging.hpp:115
std::string get_eigen_method_name(const EigenMethod &m)
Definition: naming.hpp:48
const Matrix< Scalar, Dynamic, 1 > & eigenvalues() const
Returns the eigenvalues of given matrix.
static LoggingSingleton & instance()
Definition: logging.hpp:102
EigendecompositionResult generalized_eigendecomposition_impl_arpack(const LMatrixType &lhs, const RMatrixType &rhs, IndexType target_dimension, unsigned int skip)
ARPACK implementation of eigendecomposition-based embedding.
An exception type that is thrown when unsupported method is called.
TAPKEE_INTERNAL_PAIR< tapkee::DenseMatrix, tapkee::DenseVector > EigendecompositionResult
Definition: synonyms.hpp:41
bool is(const M &m) const