Eigen  3.4.90 (git rev a4098ac676528a83cfb73d4d26ce1b42ec05f47c)
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
19#include "./InternalHeaderCheck.h"
20
21namespace Eigen {
22
23namespace internal {
24template<typename Derived> struct traits<SVDBase<Derived> >
25 : traits<Derived>
26{
27 typedef MatrixXpr XprKind;
28 typedef SolverStorage StorageKind;
29 typedef int StorageIndex;
30 enum { Flags = 0 };
31};
32}
33
64template<typename Derived> class SVDBase
65 : public SolverBase<SVDBase<Derived> >
66{
67public:
68
69 template<typename Derived_>
70 friend struct internal::solve_assertion;
71
72 typedef typename internal::traits<Derived>::MatrixType MatrixType;
73 typedef typename MatrixType::Scalar Scalar;
74 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
75 typedef typename Eigen::internal::traits<SVDBase>::StorageIndex StorageIndex;
77 enum {
78 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
79 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
80 DiagSizeAtCompileTime = internal::min_size_prefer_dynamic(RowsAtCompileTime,ColsAtCompileTime),
81 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
82 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
83 MaxDiagSizeAtCompileTime = internal::min_size_prefer_fixed(MaxRowsAtCompileTime,MaxColsAtCompileTime),
84 MatrixOptions = MatrixType::Options
85 };
86
89 typedef typename internal::plain_diag_type<MatrixType, RealScalar>::type SingularValuesType;
90
91 Derived& derived() { return *static_cast<Derived*>(this); }
92 const Derived& derived() const { return *static_cast<const Derived*>(this); }
93
103 const MatrixUType& matrixU() const
104 {
105 _check_compute_assertions();
106 eigen_assert(computeU() && "This SVD decomposition didn't compute U. Did you ask for it?");
107 return m_matrixU;
108 }
109
119 const MatrixVType& matrixV() const
120 {
121 _check_compute_assertions();
122 eigen_assert(computeV() && "This SVD decomposition didn't compute V. Did you ask for it?");
123 return m_matrixV;
124 }
125
131 const SingularValuesType& singularValues() const
132 {
133 _check_compute_assertions();
134 return m_singularValues;
135 }
136
139 {
140 _check_compute_assertions();
141 return m_nonzeroSingularValues;
142 }
143
150 inline Index rank() const
151 {
152 using std::abs;
153 _check_compute_assertions();
154 if(m_singularValues.size()==0) return 0;
155 RealScalar premultiplied_threshold = numext::maxi<RealScalar>(m_singularValues.coeff(0) * threshold(), (std::numeric_limits<RealScalar>::min)());
156 Index i = m_nonzeroSingularValues-1;
157 while(i>=0 && m_singularValues.coeff(i) < premultiplied_threshold) --i;
158 return i+1;
159 }
160
175 Derived& setThreshold(const RealScalar& threshold)
176 {
177 m_usePrescribedThreshold = true;
178 m_prescribedThreshold = threshold;
179 return derived();
180 }
181
190 Derived& setThreshold(Default_t)
191 {
192 m_usePrescribedThreshold = false;
193 return derived();
194 }
195
200 RealScalar threshold() const
201 {
202 eigen_assert(m_isInitialized || m_usePrescribedThreshold);
203 // this temporary is needed to workaround a MSVC issue
204 Index diagSize = (std::max<Index>)(1,m_diagSize);
205 return m_usePrescribedThreshold ? m_prescribedThreshold
206 : RealScalar(diagSize)*NumTraits<Scalar>::epsilon();
207 }
208
210 inline bool computeU() const { return m_computeFullU || m_computeThinU; }
212 inline bool computeV() const { return m_computeFullV || m_computeThinV; }
213
214 inline Index rows() const { return m_rows; }
215 inline Index cols() const { return m_cols; }
216
217 #ifdef EIGEN_PARSED_BY_DOXYGEN
227 template<typename Rhs>
228 inline const Solve<Derived, Rhs>
229 solve(const MatrixBase<Rhs>& b) const;
230 #endif
231
232
237 EIGEN_DEVICE_FUNC
239 {
240 eigen_assert(m_isInitialized && "SVD is not initialized.");
241 return m_info;
242 }
243
244 #ifndef EIGEN_PARSED_BY_DOXYGEN
245 template<typename RhsType, typename DstType>
246 void _solve_impl(const RhsType &rhs, DstType &dst) const;
247
248 template<bool Conjugate, typename RhsType, typename DstType>
249 void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const;
250 #endif
251
252protected:
253
254 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
255
256 void _check_compute_assertions() const {
257 eigen_assert(m_isInitialized && "SVD is not initialized.");
258 }
259
260 template<bool Transpose_, typename Rhs>
261 void _check_solve_assertion(const Rhs& b) const {
262 EIGEN_ONLY_USED_FOR_DEBUG(b);
263 _check_compute_assertions();
264 eigen_assert(computeU() && computeV() && "SVDBase::solve(): Both unitaries U and V are required to be computed (thin unitaries suffice).");
265 eigen_assert((Transpose_?cols():rows())==b.rows() && "SVDBase::solve(): invalid number of rows of the right hand side matrix b");
266 }
267
268 // return true if already allocated
269 bool allocate(Index rows, Index cols, unsigned int computationOptions) ;
270
271 MatrixUType m_matrixU;
272 MatrixVType m_matrixV;
273 SingularValuesType m_singularValues;
274 ComputationInfo m_info;
275 bool m_isInitialized, m_isAllocated, m_usePrescribedThreshold;
276 bool m_computeFullU, m_computeThinU;
277 bool m_computeFullV, m_computeThinV;
278 unsigned int m_computationOptions;
279 Index m_nonzeroSingularValues, m_rows, m_cols, m_diagSize;
280 RealScalar m_prescribedThreshold;
281
287 : m_info(Success),
288 m_isInitialized(false),
289 m_isAllocated(false),
290 m_usePrescribedThreshold(false),
291 m_computeFullU(false),
292 m_computeThinU(false),
293 m_computeFullV(false),
294 m_computeThinV(false),
295 m_computationOptions(0),
296 m_rows(-1), m_cols(-1), m_diagSize(0)
297 { }
298
299
300};
301
302#ifndef EIGEN_PARSED_BY_DOXYGEN
303template<typename Derived>
304template<typename RhsType, typename DstType>
305void SVDBase<Derived>::_solve_impl(const RhsType &rhs, DstType &dst) const
306{
307 // A = U S V^*
308 // So A^{-1} = V S^{-1} U^*
309
310 Matrix<typename RhsType::Scalar, Dynamic, RhsType::ColsAtCompileTime, 0, MatrixType::MaxRowsAtCompileTime, RhsType::MaxColsAtCompileTime> tmp;
311 Index l_rank = rank();
312 tmp.noalias() = m_matrixU.leftCols(l_rank).adjoint() * rhs;
313 tmp = m_singularValues.head(l_rank).asDiagonal().inverse() * tmp;
314 dst = m_matrixV.leftCols(l_rank) * tmp;
315}
316
317template<typename Derived>
318template<bool Conjugate, typename RhsType, typename DstType>
319void SVDBase<Derived>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
320{
321 // A = U S V^*
322 // So A^{-*} = U S^{-1} V^*
323 // And A^{-T} = U_conj S^{-1} V^T
324 Matrix<typename RhsType::Scalar, Dynamic, RhsType::ColsAtCompileTime, 0, MatrixType::MaxRowsAtCompileTime, RhsType::MaxColsAtCompileTime> tmp;
325 Index l_rank = rank();
326
327 tmp.noalias() = m_matrixV.leftCols(l_rank).transpose().template conjugateIf<Conjugate>() * rhs;
328 tmp = m_singularValues.head(l_rank).asDiagonal().inverse() * tmp;
329 dst = m_matrixU.template conjugateIf<!Conjugate>().leftCols(l_rank) * tmp;
330}
331#endif
332
333template<typename MatrixType>
334bool SVDBase<MatrixType>::allocate(Index rows, Index cols, unsigned int computationOptions)
335{
336 eigen_assert(rows >= 0 && cols >= 0);
337
338 if (m_isAllocated &&
339 rows == m_rows &&
340 cols == m_cols &&
341 computationOptions == m_computationOptions)
342 {
343 return true;
344 }
345
346 m_rows = rows;
347 m_cols = cols;
348 m_info = Success;
349 m_isInitialized = false;
350 m_isAllocated = true;
351 m_computationOptions = computationOptions;
352 m_computeFullU = (computationOptions & ComputeFullU) != 0;
353 m_computeThinU = (computationOptions & ComputeThinU) != 0;
354 m_computeFullV = (computationOptions & ComputeFullV) != 0;
355 m_computeThinV = (computationOptions & ComputeThinV) != 0;
356 eigen_assert(!(m_computeFullU && m_computeThinU) && "SVDBase: you can't ask for both full and thin U");
357 eigen_assert(!(m_computeFullV && m_computeThinV) && "SVDBase: you can't ask for both full and thin V");
358 eigen_assert(internal::check_implication(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==Dynamic) &&
359 "SVDBase: thin U and V are only available when your matrix has a dynamic number of columns.");
360
361 m_diagSize = (std::min)(m_rows, m_cols);
362 m_singularValues.resize(m_diagSize);
363 if(RowsAtCompileTime==Dynamic)
364 m_matrixU.resize(m_rows, m_computeFullU ? m_rows : m_computeThinU ? m_diagSize : 0);
365 if(ColsAtCompileTime==Dynamic)
366 m_matrixV.resize(m_cols, m_computeFullV ? m_cols : m_computeThinV ? m_diagSize : 0);
367
368 return false;
369}
370
371}// end namespace
372
373#endif // EIGEN_SVDBASE_H
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:52
Base class of SVD algorithms.
Definition: SVDBase.h:66
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SVDBase.h:238
Derived & setThreshold(const RealScalar &threshold)
Definition: SVDBase.h:175
const Solve< Derived, Rhs > solve(const MatrixBase< Rhs > &b) const
Index rank() const
Definition: SVDBase.h:150
bool computeV() const
Definition: SVDBase.h:212
Eigen::Index Index
Definition: SVDBase.h:76
bool computeU() const
Definition: SVDBase.h:210
Derived & setThreshold(Default_t)
Definition: SVDBase.h:190
RealScalar threshold() const
Definition: SVDBase.h:200
SVDBase()
Default Constructor.
Definition: SVDBase.h:286
const SingularValuesType & singularValues() const
Definition: SVDBase.h:131
const MatrixUType & matrixU() const
Definition: SVDBase.h:103
const MatrixVType & matrixV() const
Definition: SVDBase.h:119
Index nonzeroSingularValues() const
Definition: SVDBase.h:138
A base class for matrix decomposition and solvers.
Definition: SolverBase.h:71
ComputationInfo
Definition: Constants.h:442
@ Success
Definition: Constants.h:444
@ ComputeFullV
Definition: Constants.h:399
@ ComputeThinV
Definition: Constants.h:401
@ ComputeFullU
Definition: Constants.h:395
@ ComputeThinU
Definition: Constants.h:397
Namespace containing all symbols from the Eigen library.
Definition: B01_Experimental.dox:1
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
const int Dynamic
Definition: Constants.h:24
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:41
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:235