Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
Eigen  3.4.0
 
Loading...
Searching...
No Matches
SVDBase.h
1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2009-2010 Benoit Jacob <jacob.benoit.1@gmail.com>
5// Copyright (C) 2014 Gael Guennebaud <gael.guennebaud@inria.fr>
6//
7// Copyright (C) 2013 Gauthier Brun <brun.gauthier@gmail.com>
8// Copyright (C) 2013 Nicolas Carre <nicolas.carre@ensimag.fr>
9// Copyright (C) 2013 Jean Ceccato <jean.ceccato@ensimag.fr>
10// Copyright (C) 2013 Pierre Zoppitelli <pierre.zoppitelli@ensimag.fr>
11//
12// This Source Code Form is subject to the terms of the Mozilla
13// Public License v. 2.0. If a copy of the MPL was not distributed
14// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
15
16#ifndef EIGEN_SVDBASE_H
17#define EIGEN_SVDBASE_H
18
19namespace Eigen {
20
21namespace internal {
22template<typename Derived> struct traits<SVDBase<Derived> >
23 : traits<Derived>
24{
25 typedef MatrixXpr XprKind;
26 typedef SolverStorage StorageKind;
27 typedef int StorageIndex;
28 enum { Flags = 0 };
29};
30}
31
62template<typename Derived> class SVDBase
63 : public SolverBase<SVDBase<Derived> >
64{
65public:
66
67 template<typename Derived_>
68 friend struct internal::solve_assertion;
69
70 typedef typename internal::traits<Derived>::MatrixType MatrixType;
71 typedef typename MatrixType::Scalar Scalar;
72 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
73 typedef typename Eigen::internal::traits<SVDBase>::StorageIndex StorageIndex;
75 enum {
76 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
77 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
78 DiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime,ColsAtCompileTime),
79 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
80 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
81 MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(MaxRowsAtCompileTime,MaxColsAtCompileTime),
82 MatrixOptions = MatrixType::Options
83 };
84
87 typedef typename internal::plain_diag_type<MatrixType, RealScalar>::type SingularValuesType;
88
89 Derived& derived() { return *static_cast<Derived*>(this); }
90 const Derived& derived() const { return *static_cast<const Derived*>(this); }
91
101 const MatrixUType& matrixU() const
102 {
103 _check_compute_assertions();
104 eigen_assert(computeU() && "This SVD decomposition didn't compute U. Did you ask for it?");
105 return m_matrixU;
106 }
107
117 const MatrixVType& matrixV() const
118 {
119 _check_compute_assertions();
120 eigen_assert(computeV() && "This SVD decomposition didn't compute V. Did you ask for it?");
121 return m_matrixV;
122 }
123
129 const SingularValuesType& singularValues() const
130 {
131 _check_compute_assertions();
132 return m_singularValues;
133 }
134
137 {
138 _check_compute_assertions();
139 return m_nonzeroSingularValues;
140 }
141
148 inline Index rank() const
149 {
150 using std::abs;
151 _check_compute_assertions();
152 if(m_singularValues.size()==0) return 0;
153 RealScalar premultiplied_threshold = numext::maxi<RealScalar>(m_singularValues.coeff(0) * threshold(), (std::numeric_limits<RealScalar>::min)());
154 Index i = m_nonzeroSingularValues-1;
155 while(i>=0 && m_singularValues.coeff(i) < premultiplied_threshold) --i;
156 return i+1;
157 }
158
173 Derived& setThreshold(const RealScalar& threshold)
174 {
175 m_usePrescribedThreshold = true;
176 m_prescribedThreshold = threshold;
177 return derived();
178 }
179
188 Derived& setThreshold(Default_t)
189 {
190 m_usePrescribedThreshold = false;
191 return derived();
192 }
193
198 RealScalar threshold() const
199 {
200 eigen_assert(m_isInitialized || m_usePrescribedThreshold);
201 // this temporary is needed to workaround a MSVC issue
202 Index diagSize = (std::max<Index>)(1,m_diagSize);
203 return m_usePrescribedThreshold ? m_prescribedThreshold
204 : RealScalar(diagSize)*NumTraits<Scalar>::epsilon();
205 }
206
208 inline bool computeU() const { return m_computeFullU || m_computeThinU; }
210 inline bool computeV() const { return m_computeFullV || m_computeThinV; }
211
212 inline Index rows() const { return m_rows; }
213 inline Index cols() const { return m_cols; }
214
215 #ifdef EIGEN_PARSED_BY_DOXYGEN
225 template<typename Rhs>
226 inline const Solve<Derived, Rhs>
227 solve(const MatrixBase<Rhs>& b) const;
228 #endif
229
230
235 EIGEN_DEVICE_FUNC
237 {
238 eigen_assert(m_isInitialized && "SVD is not initialized.");
239 return m_info;
240 }
241
242 #ifndef EIGEN_PARSED_BY_DOXYGEN
243 template<typename RhsType, typename DstType>
244 void _solve_impl(const RhsType &rhs, DstType &dst) const;
245
246 template<bool Conjugate, typename RhsType, typename DstType>
247 void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const;
248 #endif
249
250protected:
251
252 static void check_template_parameters()
253 {
254 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
255 }
256
257 void _check_compute_assertions() const {
258 eigen_assert(m_isInitialized && "SVD is not initialized.");
259 }
260
261 template<bool Transpose_, typename Rhs>
262 void _check_solve_assertion(const Rhs& b) const {
263 EIGEN_ONLY_USED_FOR_DEBUG(b);
264 _check_compute_assertions();
265 eigen_assert(computeU() && computeV() && "SVDBase::solve(): Both unitaries U and V are required to be computed (thin unitaries suffice).");
266 eigen_assert((Transpose_?cols():rows())==b.rows() && "SVDBase::solve(): invalid number of rows of the right hand side matrix b");
267 }
268
269 // return true if already allocated
270 bool allocate(Index rows, Index cols, unsigned int computationOptions) ;
271
272 MatrixUType m_matrixU;
273 MatrixVType m_matrixV;
274 SingularValuesType m_singularValues;
275 ComputationInfo m_info;
276 bool m_isInitialized, m_isAllocated, m_usePrescribedThreshold;
277 bool m_computeFullU, m_computeThinU;
278 bool m_computeFullV, m_computeThinV;
279 unsigned int m_computationOptions;
280 Index m_nonzeroSingularValues, m_rows, m_cols, m_diagSize;
281 RealScalar m_prescribedThreshold;
282
288 : m_info(Success),
289 m_isInitialized(false),
290 m_isAllocated(false),
291 m_usePrescribedThreshold(false),
292 m_computeFullU(false),
293 m_computeThinU(false),
294 m_computeFullV(false),
295 m_computeThinV(false),
296 m_computationOptions(0),
297 m_rows(-1), m_cols(-1), m_diagSize(0)
298 {
299 check_template_parameters();
300 }
301
302
303};
304
305#ifndef EIGEN_PARSED_BY_DOXYGEN
306template<typename Derived>
307template<typename RhsType, typename DstType>
308void SVDBase<Derived>::_solve_impl(const RhsType &rhs, DstType &dst) const
309{
310 // A = U S V^*
311 // So A^{-1} = V S^{-1} U^*
312
313 Matrix<typename RhsType::Scalar, Dynamic, RhsType::ColsAtCompileTime, 0, MatrixType::MaxRowsAtCompileTime, RhsType::MaxColsAtCompileTime> tmp;
314 Index l_rank = rank();
315 tmp.noalias() = m_matrixU.leftCols(l_rank).adjoint() * rhs;
316 tmp = m_singularValues.head(l_rank).asDiagonal().inverse() * tmp;
317 dst = m_matrixV.leftCols(l_rank) * tmp;
318}
319
320template<typename Derived>
321template<bool Conjugate, typename RhsType, typename DstType>
322void SVDBase<Derived>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
323{
324 // A = U S V^*
325 // So A^{-*} = U S^{-1} V^*
326 // And A^{-T} = U_conj S^{-1} V^T
327 Matrix<typename RhsType::Scalar, Dynamic, RhsType::ColsAtCompileTime, 0, MatrixType::MaxRowsAtCompileTime, RhsType::MaxColsAtCompileTime> tmp;
328 Index l_rank = rank();
329
330 tmp.noalias() = m_matrixV.leftCols(l_rank).transpose().template conjugateIf<Conjugate>() * rhs;
331 tmp = m_singularValues.head(l_rank).asDiagonal().inverse() * tmp;
332 dst = m_matrixU.template conjugateIf<!Conjugate>().leftCols(l_rank) * tmp;
333}
334#endif
335
336template<typename MatrixType>
337bool SVDBase<MatrixType>::allocate(Index rows, Index cols, unsigned int computationOptions)
338{
339 eigen_assert(rows >= 0 && cols >= 0);
340
341 if (m_isAllocated &&
342 rows == m_rows &&
343 cols == m_cols &&
344 computationOptions == m_computationOptions)
345 {
346 return true;
347 }
348
349 m_rows = rows;
350 m_cols = cols;
351 m_info = Success;
352 m_isInitialized = false;
353 m_isAllocated = true;
354 m_computationOptions = computationOptions;
355 m_computeFullU = (computationOptions & ComputeFullU) != 0;
356 m_computeThinU = (computationOptions & ComputeThinU) != 0;
357 m_computeFullV = (computationOptions & ComputeFullV) != 0;
358 m_computeThinV = (computationOptions & ComputeThinV) != 0;
359 eigen_assert(!(m_computeFullU && m_computeThinU) && "SVDBase: you can't ask for both full and thin U");
360 eigen_assert(!(m_computeFullV && m_computeThinV) && "SVDBase: you can't ask for both full and thin V");
361 eigen_assert(EIGEN_IMPLIES(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==Dynamic) &&
362 "SVDBase: thin U and V are only available when your matrix has a dynamic number of columns.");
363
364 m_diagSize = (std::min)(m_rows, m_cols);
365 m_singularValues.resize(m_diagSize);
366 if(RowsAtCompileTime==Dynamic)
367 m_matrixU.resize(m_rows, m_computeFullU ? m_rows : m_computeThinU ? m_diagSize : 0);
368 if(ColsAtCompileTime==Dynamic)
369 m_matrixV.resize(m_cols, m_computeFullV ? m_cols : m_computeThinV ? m_diagSize : 0);
370
371 return false;
372}
373
374}// end namespace
375
376#endif // EIGEN_SVDBASE_H
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:50
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:180
Base class of SVD algorithms.
Definition: SVDBase.h:64
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SVDBase.h:236
Derived & setThreshold(const RealScalar &threshold)
Definition: SVDBase.h:173
const Solve< Derived, Rhs > solve(const MatrixBase< Rhs > &b) const
Index rank() const
Definition: SVDBase.h:148
bool computeV() const
Definition: SVDBase.h:210
Eigen::Index Index
Definition: SVDBase.h:74
bool computeU() const
Definition: SVDBase.h:208
Derived & setThreshold(Default_t)
Definition: SVDBase.h:188
RealScalar threshold() const
Definition: SVDBase.h:198
SVDBase()
Default Constructor.
Definition: SVDBase.h:287
const SingularValuesType & singularValues() const
Definition: SVDBase.h:129
const MatrixUType & matrixU() const
Definition: SVDBase.h:101
const MatrixVType & matrixV() const
Definition: SVDBase.h:117
Index nonzeroSingularValues() const
Definition: SVDBase.h:136
A base class for matrix decomposition and solvers.
Definition: SolverBase.h:69
ComputationInfo
Definition: Constants.h:440
@ Success
Definition: Constants.h:442
@ ComputeFullV
Definition: Constants.h:397
@ ComputeThinV
Definition: Constants.h:399
@ ComputeFullU
Definition: Constants.h:393
@ ComputeThinU
Definition: Constants.h:395
Namespace containing all symbols from the Eigen library.
Definition: Core:141
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:233