Eigen  3.3.3
SVDBase.h
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2009-2010 Benoit Jacob <jacob.benoit.1@gmail.com>
00005 // Copyright (C) 2014 Gael Guennebaud <gael.guennebaud@inria.fr>
00006 //
00007 // Copyright (C) 2013 Gauthier Brun <brun.gauthier@gmail.com>
00008 // Copyright (C) 2013 Nicolas Carre <nicolas.carre@ensimag.fr>
00009 // Copyright (C) 2013 Jean Ceccato <jean.ceccato@ensimag.fr>
00010 // Copyright (C) 2013 Pierre Zoppitelli <pierre.zoppitelli@ensimag.fr>
00011 //
00012 // This Source Code Form is subject to the terms of the Mozilla
00013 // Public License v. 2.0. If a copy of the MPL was not distributed
00014 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00015 
00016 #ifndef EIGEN_SVDBASE_H
00017 #define EIGEN_SVDBASE_H
00018 
00019 namespace Eigen {
00047 template<typename Derived>
00048 class SVDBase
00049 {
00050 
00051 public:
00052   typedef typename internal::traits<Derived>::MatrixType MatrixType;
00053   typedef typename MatrixType::Scalar Scalar;
00054   typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
00055   typedef typename MatrixType::StorageIndex StorageIndex;
00056   typedef Eigen::Index Index; 
00057   enum {
00058     RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00059     ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00060     DiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime,ColsAtCompileTime),
00061     MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00062     MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
00063     MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(MaxRowsAtCompileTime,MaxColsAtCompileTime),
00064     MatrixOptions = MatrixType::Options
00065   };
00066 
00067   typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, MatrixOptions, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixUType;
00068   typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime, MatrixOptions, MaxColsAtCompileTime, MaxColsAtCompileTime> MatrixVType;
00069   typedef typename internal::plain_diag_type<MatrixType, RealScalar>::type SingularValuesType;
00070   
00071   Derived& derived() { return *static_cast<Derived*>(this); }
00072   const Derived& derived() const { return *static_cast<const Derived*>(this); }
00073 
00083   const MatrixUType& matrixU() const
00084   {
00085     eigen_assert(m_isInitialized && "SVD is not initialized.");
00086     eigen_assert(computeU() && "This SVD decomposition didn't compute U. Did you ask for it?");
00087     return m_matrixU;
00088   }
00089 
00099   const MatrixVType& matrixV() const
00100   {
00101     eigen_assert(m_isInitialized && "SVD is not initialized.");
00102     eigen_assert(computeV() && "This SVD decomposition didn't compute V. Did you ask for it?");
00103     return m_matrixV;
00104   }
00105 
00111   const SingularValuesType& singularValues() const
00112   {
00113     eigen_assert(m_isInitialized && "SVD is not initialized.");
00114     return m_singularValues;
00115   }
00116 
00118   Index nonzeroSingularValues() const
00119   {
00120     eigen_assert(m_isInitialized && "SVD is not initialized.");
00121     return m_nonzeroSingularValues;
00122   }
00123   
00130   inline Index rank() const
00131   {
00132     using std::abs;
00133     eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
00134     if(m_singularValues.size()==0) return 0;
00135     RealScalar premultiplied_threshold = numext::maxi<RealScalar>(m_singularValues.coeff(0) * threshold(), (std::numeric_limits<RealScalar>::min)());
00136     Index i = m_nonzeroSingularValues-1;
00137     while(i>=0 && m_singularValues.coeff(i) < premultiplied_threshold) --i;
00138     return i+1;
00139   }
00140   
00155   Derived& setThreshold(const RealScalar& threshold)
00156   {
00157     m_usePrescribedThreshold = true;
00158     m_prescribedThreshold = threshold;
00159     return derived();
00160   }
00161 
00170   Derived& setThreshold(Default_t)
00171   {
00172     m_usePrescribedThreshold = false;
00173     return derived();
00174   }
00175 
00180   RealScalar threshold() const
00181   {
00182     eigen_assert(m_isInitialized || m_usePrescribedThreshold);
00183     return m_usePrescribedThreshold ? m_prescribedThreshold
00184                                     : (std::max<Index>)(1,m_diagSize)*NumTraits<Scalar>::epsilon();
00185   }
00186 
00188   inline bool computeU() const { return m_computeFullU || m_computeThinU; }
00190   inline bool computeV() const { return m_computeFullV || m_computeThinV; }
00191 
00192   inline Index rows() const { return m_rows; }
00193   inline Index cols() const { return m_cols; }
00194   
00204   template<typename Rhs>
00205   inline const Solve<Derived, Rhs>
00206   solve(const MatrixBase<Rhs>& b) const
00207   {
00208     eigen_assert(m_isInitialized && "SVD is not initialized.");
00209     eigen_assert(computeU() && computeV() && "SVD::solve() requires both unitaries U and V to be computed (thin unitaries suffice).");
00210     return Solve<Derived, Rhs>(derived(), b.derived());
00211   }
00212   
00213   #ifndef EIGEN_PARSED_BY_DOXYGEN
00214   template<typename RhsType, typename DstType>
00215   EIGEN_DEVICE_FUNC
00216   void _solve_impl(const RhsType &rhs, DstType &dst) const;
00217   #endif
00218 
00219 protected:
00220   
00221   static void check_template_parameters()
00222   {
00223     EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
00224   }
00225   
00226   // return true if already allocated
00227   bool allocate(Index rows, Index cols, unsigned int computationOptions) ;
00228 
00229   MatrixUType m_matrixU;
00230   MatrixVType m_matrixV;
00231   SingularValuesType m_singularValues;
00232   bool m_isInitialized, m_isAllocated, m_usePrescribedThreshold;
00233   bool m_computeFullU, m_computeThinU;
00234   bool m_computeFullV, m_computeThinV;
00235   unsigned int m_computationOptions;
00236   Index m_nonzeroSingularValues, m_rows, m_cols, m_diagSize;
00237   RealScalar m_prescribedThreshold;
00238 
00243   SVDBase()
00244     : m_isInitialized(false),
00245       m_isAllocated(false),
00246       m_usePrescribedThreshold(false),
00247       m_computationOptions(0),
00248       m_rows(-1), m_cols(-1), m_diagSize(0)
00249   {
00250     check_template_parameters();
00251   }
00252 
00253 
00254 };
00255 
00256 #ifndef EIGEN_PARSED_BY_DOXYGEN
00257 template<typename Derived>
00258 template<typename RhsType, typename DstType>
00259 void SVDBase<Derived>::_solve_impl(const RhsType &rhs, DstType &dst) const
00260 {
00261   eigen_assert(rhs.rows() == rows());
00262 
00263   // A = U S V^*
00264   // So A^{-1} = V S^{-1} U^*
00265 
00266   Matrix<Scalar, Dynamic, RhsType::ColsAtCompileTime, 0, MatrixType::MaxRowsAtCompileTime, RhsType::MaxColsAtCompileTime> tmp;
00267   Index l_rank = rank();
00268   tmp.noalias() =  m_matrixU.leftCols(l_rank).adjoint() * rhs;
00269   tmp = m_singularValues.head(l_rank).asDiagonal().inverse() * tmp;
00270   dst = m_matrixV.leftCols(l_rank) * tmp;
00271 }
00272 #endif
00273 
00274 template<typename MatrixType>
00275 bool SVDBase<MatrixType>::allocate(Index rows, Index cols, unsigned int computationOptions)
00276 {
00277   eigen_assert(rows >= 0 && cols >= 0);
00278 
00279   if (m_isAllocated &&
00280       rows == m_rows &&
00281       cols == m_cols &&
00282       computationOptions == m_computationOptions)
00283   {
00284     return true;
00285   }
00286 
00287   m_rows = rows;
00288   m_cols = cols;
00289   m_isInitialized = false;
00290   m_isAllocated = true;
00291   m_computationOptions = computationOptions;
00292   m_computeFullU = (computationOptions & ComputeFullU) != 0;
00293   m_computeThinU = (computationOptions & ComputeThinU) != 0;
00294   m_computeFullV = (computationOptions & ComputeFullV) != 0;
00295   m_computeThinV = (computationOptions & ComputeThinV) != 0;
00296   eigen_assert(!(m_computeFullU && m_computeThinU) && "SVDBase: you can't ask for both full and thin U");
00297   eigen_assert(!(m_computeFullV && m_computeThinV) && "SVDBase: you can't ask for both full and thin V");
00298   eigen_assert(EIGEN_IMPLIES(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==Dynamic) &&
00299                "SVDBase: thin U and V are only available when your matrix has a dynamic number of columns.");
00300 
00301   m_diagSize = (std::min)(m_rows, m_cols);
00302   m_singularValues.resize(m_diagSize);
00303   if(RowsAtCompileTime==Dynamic)
00304     m_matrixU.resize(m_rows, m_computeFullU ? m_rows : m_computeThinU ? m_diagSize : 0);
00305   if(ColsAtCompileTime==Dynamic)
00306     m_matrixV.resize(m_cols, m_computeFullV ? m_cols : m_computeThinV ? m_diagSize : 0);
00307 
00308   return false;
00309 }
00310 
00311 }// end namespace
00312 
00313 #endif // EIGEN_SVDBASE_H
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends