diff options
Diffstat (limited to 'Eigen/src/SVD/JacobiSVD.h')
| -rw-r--r-- | Eigen/src/SVD/JacobiSVD.h | 246 |
1 files changed, 126 insertions, 120 deletions
diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h index 56ffea131..d180db7db 100644 --- a/Eigen/src/SVD/JacobiSVD.h +++ b/Eigen/src/SVD/JacobiSVD.h @@ -25,12 +25,12 @@ #ifndef EIGEN_JACOBISVD_H #define EIGEN_JACOBISVD_H +namespace internal { // forward declaration (needed by ICC) // the empty body is required by MSVC template<typename MatrixType, int QRPreconditioner, bool IsComplex = NumTraits<typename MatrixType::Scalar>::IsComplex> -struct ei_svd_precondition_2x2_block_to_be_real {}; - +struct svd_precondition_2x2_block_to_be_real {}; /*** QR preconditioners (R-SVD) *** @@ -42,7 +42,7 @@ struct ei_svd_precondition_2x2_block_to_be_real {}; enum { PreconditionIfMoreColsThanRows, PreconditionIfMoreRowsThanCols }; template<typename MatrixType, int QRPreconditioner, int Case> -struct ei_qr_preconditioner_should_do_anything +struct qr_preconditioner_should_do_anything { enum { a = MatrixType::RowsAtCompileTime != Dynamic && MatrixType::ColsAtCompileTime != Dynamic && @@ -57,11 +57,11 @@ struct ei_qr_preconditioner_should_do_anything }; template<typename MatrixType, int QRPreconditioner, int Case, - bool DoAnything = ei_qr_preconditioner_should_do_anything<MatrixType, QRPreconditioner, Case>::ret -> struct ei_qr_preconditioner_impl {}; + bool DoAnything = qr_preconditioner_should_do_anything<MatrixType, QRPreconditioner, Case>::ret +> struct qr_preconditioner_impl {}; template<typename MatrixType, int QRPreconditioner, int Case> -struct ei_qr_preconditioner_impl<MatrixType, QRPreconditioner, Case, false> +struct qr_preconditioner_impl<MatrixType, QRPreconditioner, Case, false> { static bool run(JacobiSVD<MatrixType, QRPreconditioner>&, const MatrixType&) { @@ -72,7 +72,7 @@ struct ei_qr_preconditioner_impl<MatrixType, QRPreconditioner, Case, false> /*** preconditioner using FullPivHouseholderQR ***/ template<typename MatrixType> -struct ei_qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true> +struct qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true> { static bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix) { @@ -89,7 +89,7 @@ struct ei_qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, }; template<typename MatrixType> -struct ei_qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true> +struct qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true> { static bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix) { @@ -111,7 +111,7 @@ struct ei_qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, /*** preconditioner using ColPivHouseholderQR ***/ template<typename MatrixType> -struct ei_qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true> +struct qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true> { static bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix) { @@ -132,7 +132,7 @@ struct ei_qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, }; template<typename MatrixType> -struct ei_qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true> +struct qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true> { static bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix) { @@ -158,7 +158,7 @@ struct ei_qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, /*** preconditioner using HouseholderQR ***/ template<typename MatrixType> -struct ei_qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true> +struct qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true> { static bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix) { @@ -179,7 +179,7 @@ struct ei_qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, Precon }; template<typename MatrixType> -struct ei_qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true> +struct qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true> { static bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix) { @@ -202,7 +202,90 @@ struct ei_qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, Precon } }; +/*** 2x2 SVD implementation + *** + *** JacobiSVD consists in performing a series of 2x2 SVD subproblems + ***/ + +template<typename MatrixType, int QRPreconditioner> +struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, false> +{ + typedef JacobiSVD<MatrixType, QRPreconditioner> SVD; + typedef typename SVD::Index Index; + static void run(typename SVD::WorkMatrixType&, SVD&, Index, Index) {} +}; + +template<typename MatrixType, int QRPreconditioner> +struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true> +{ + typedef JacobiSVD<MatrixType, QRPreconditioner> SVD; + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef typename SVD::Index Index; + static void run(typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index p, Index q) + { + Scalar z; + JacobiRotation<Scalar> rot; + RealScalar n = sqrt(abs2(work_matrix.coeff(p,p)) + abs2(work_matrix.coeff(q,p))); + if(n==0) + { + z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q); + work_matrix.row(p) *= z; + if(svd.computeU()) svd.m_matrixU.col(p) *= conj(z); + z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q); + work_matrix.row(q) *= z; + if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z); + } + else + { + rot.c() = conj(work_matrix.coeff(p,p)) / n; + rot.s() = work_matrix.coeff(q,p) / n; + work_matrix.applyOnTheLeft(p,q,rot); + if(svd.computeU()) svd.m_matrixU.applyOnTheRight(p,q,rot.adjoint()); + if(work_matrix.coeff(p,q) != Scalar(0)) + { + Scalar z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q); + work_matrix.col(q) *= z; + if(svd.computeV()) svd.m_matrixV.col(q) *= z; + } + if(work_matrix.coeff(q,q) != Scalar(0)) + { + z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q); + work_matrix.row(q) *= z; + if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z); + } + } + } +}; + +template<typename MatrixType, typename RealScalar, typename Index> +void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q, + JacobiRotation<RealScalar> *j_left, + JacobiRotation<RealScalar> *j_right) +{ + Matrix<RealScalar,2,2> m; + m << real(matrix.coeff(p,p)), real(matrix.coeff(p,q)), + real(matrix.coeff(q,p)), real(matrix.coeff(q,q)); + JacobiRotation<RealScalar> rot1; + RealScalar t = m.coeff(0,0) + m.coeff(1,1); + RealScalar d = m.coeff(1,0) - m.coeff(0,1); + if(t == RealScalar(0)) + { + rot1.c() = 0; + rot1.s() = d > 0 ? 1 : -1; + } + else + { + RealScalar u = d / t; + rot1.c() = RealScalar(1) / sqrt(1 + abs2(u)); + rot1.s() = rot1.c() * u; + } + m.applyOnTheLeft(0,1,rot1); + j_right->makeJacobi(m,0,1); + *j_left = rot1 * j_right->transpose(); +} +} // end namespace internal /** \ingroup SVD_Module * @@ -281,9 +364,9 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime, MatrixOptions, MaxColsAtCompileTime, MaxColsAtCompileTime> MatrixVType; - typedef typename ei_plain_diag_type<MatrixType, RealScalar>::type SingularValuesType; - typedef typename ei_plain_row_type<MatrixType>::type RowType; - typedef typename ei_plain_col_type<MatrixType>::type ColType; + typedef typename internal::plain_diag_type<MatrixType, RealScalar>::type SingularValuesType; + typedef typename internal::plain_row_type<MatrixType>::type RowType; + typedef typename internal::plain_col_type<MatrixType>::type ColType; typedef Matrix<Scalar, DiagSizeAtCompileTime, DiagSizeAtCompileTime, MatrixOptions, MaxDiagSizeAtCompileTime, MaxDiagSizeAtCompileTime> WorkMatrixType; @@ -345,8 +428,8 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD */ const MatrixUType& matrixU() const { - ei_assert(m_isInitialized && "JacobiSVD is not initialized."); - ei_assert(computeU() && "This JacobiSVD decomposition didn't compute U. Did you ask for it?"); + eigen_assert(m_isInitialized && "JacobiSVD is not initialized."); + eigen_assert(computeU() && "This JacobiSVD decomposition didn't compute U. Did you ask for it?"); return m_matrixU; } @@ -361,8 +444,8 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD */ const MatrixVType& matrixV() const { - ei_assert(m_isInitialized && "JacobiSVD is not initialized."); - ei_assert(computeV() && "This JacobiSVD decomposition didn't compute V. Did you ask for it?"); + eigen_assert(m_isInitialized && "JacobiSVD is not initialized."); + eigen_assert(computeV() && "This JacobiSVD decomposition didn't compute V. Did you ask for it?"); return m_matrixV; } @@ -373,7 +456,7 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD */ const SingularValuesType& singularValues() const { - ei_assert(m_isInitialized && "JacobiSVD is not initialized."); + eigen_assert(m_isInitialized && "JacobiSVD is not initialized."); return m_singularValues; } @@ -392,18 +475,18 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD * In other words, the returned solution is guaranteed to minimize the Euclidean norm \f$ \Vert A x - b \Vert \f$. */ template<typename Rhs> - inline const ei_solve_retval<JacobiSVD, Rhs> + inline const internal::solve_retval<JacobiSVD, Rhs> solve(const MatrixBase<Rhs>& b) const { - ei_assert(m_isInitialized && "JacobiSVD is not initialized."); - ei_assert(computeU() && computeV() && "JacobiSVD::solve() requires both unitaries U and V to be computed (thin unitaries suffice)."); - return ei_solve_retval<JacobiSVD, Rhs>(*this, b.derived()); + eigen_assert(m_isInitialized && "JacobiSVD is not initialized."); + eigen_assert(computeU() && computeV() && "JacobiSVD::solve() requires both unitaries U and V to be computed (thin unitaries suffice)."); + return internal::solve_retval<JacobiSVD, Rhs>(*this, b.derived()); } /** \returns the number of singular values that are not exactly 0 */ Index nonzeroSingularValues() const { - ei_assert(m_isInitialized && "JacobiSVD is not initialized."); + eigen_assert(m_isInitialized && "JacobiSVD is not initialized."); return m_nonzeroSingularValues; } @@ -424,9 +507,9 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD Index m_nonzeroSingularValues, m_rows, m_cols, m_diagSize; template<typename __MatrixType, int _QRPreconditioner, bool _IsComplex> - friend struct ei_svd_precondition_2x2_block_to_be_real; + friend struct internal::svd_precondition_2x2_block_to_be_real; template<typename __MatrixType, int _QRPreconditioner, int _Case, bool _DoAnything> - friend struct ei_qr_preconditioner_impl; + friend struct internal::qr_preconditioner_impl; }; template<typename MatrixType, int QRPreconditioner> @@ -439,13 +522,13 @@ void JacobiSVD<MatrixType, QRPreconditioner>::allocate(Index rows, Index cols, u m_computeThinU = (computationOptions & ComputeThinU) != 0; m_computeFullV = (computationOptions & ComputeFullV) != 0; m_computeThinV = (computationOptions & ComputeThinV) != 0; - ei_assert(!(m_computeFullU && m_computeThinU) && "JacobiSVD: you can't ask for both full and thin U"); - ei_assert(!(m_computeFullV && m_computeThinV) && "JacobiSVD: you can't ask for both full and thin V"); - ei_assert(EIGEN_IMPLIES(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==Dynamic) && + eigen_assert(!(m_computeFullU && m_computeThinU) && "JacobiSVD: you can't ask for both full and thin U"); + eigen_assert(!(m_computeFullV && m_computeThinV) && "JacobiSVD: you can't ask for both full and thin V"); + eigen_assert(EIGEN_IMPLIES(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==Dynamic) && "JacobiSVD: thin U and V are only available when your matrix has a dynamic number of columns."); if (QRPreconditioner == FullPivHouseholderQRPreconditioner) { - ei_assert(!(m_computeThinU || m_computeThinV) && + eigen_assert(!(m_computeThinU || m_computeThinV) && "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. " "Use the ColPivHouseholderQR preconditioner instead."); } @@ -460,85 +543,6 @@ void JacobiSVD<MatrixType, QRPreconditioner>::allocate(Index rows, Index cols, u m_workMatrix.resize(m_diagSize, m_diagSize); } - -template<typename MatrixType, int QRPreconditioner> -struct ei_svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, false> -{ - typedef JacobiSVD<MatrixType, QRPreconditioner> SVD; - typedef typename SVD::Index Index; - static void run(typename SVD::WorkMatrixType&, SVD&, Index, Index) {} -}; - -template<typename MatrixType, int QRPreconditioner> -struct ei_svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true> -{ - typedef JacobiSVD<MatrixType, QRPreconditioner> SVD; - typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::RealScalar RealScalar; - typedef typename SVD::Index Index; - static void run(typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index p, Index q) - { - Scalar z; - JacobiRotation<Scalar> rot; - RealScalar n = ei_sqrt(ei_abs2(work_matrix.coeff(p,p)) + ei_abs2(work_matrix.coeff(q,p))); - if(n==0) - { - z = ei_abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q); - work_matrix.row(p) *= z; - if(svd.computeU()) svd.m_matrixU.col(p) *= ei_conj(z); - z = ei_abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q); - work_matrix.row(q) *= z; - if(svd.computeU()) svd.m_matrixU.col(q) *= ei_conj(z); - } - else - { - rot.c() = ei_conj(work_matrix.coeff(p,p)) / n; - rot.s() = work_matrix.coeff(q,p) / n; - work_matrix.applyOnTheLeft(p,q,rot); - if(svd.computeU()) svd.m_matrixU.applyOnTheRight(p,q,rot.adjoint()); - if(work_matrix.coeff(p,q) != Scalar(0)) - { - Scalar z = ei_abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q); - work_matrix.col(q) *= z; - if(svd.computeV()) svd.m_matrixV.col(q) *= z; - } - if(work_matrix.coeff(q,q) != Scalar(0)) - { - z = ei_abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q); - work_matrix.row(q) *= z; - if(svd.computeU()) svd.m_matrixU.col(q) *= ei_conj(z); - } - } - } -}; - -template<typename MatrixType, typename RealScalar, typename Index> -void ei_real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q, - JacobiRotation<RealScalar> *j_left, - JacobiRotation<RealScalar> *j_right) -{ - Matrix<RealScalar,2,2> m; - m << ei_real(matrix.coeff(p,p)), ei_real(matrix.coeff(p,q)), - ei_real(matrix.coeff(q,p)), ei_real(matrix.coeff(q,q)); - JacobiRotation<RealScalar> rot1; - RealScalar t = m.coeff(0,0) + m.coeff(1,1); - RealScalar d = m.coeff(1,0) - m.coeff(0,1); - if(t == RealScalar(0)) - { - rot1.c() = 0; - rot1.s() = d > 0 ? 1 : -1; - } - else - { - RealScalar u = d / t; - rot1.c() = RealScalar(1) / ei_sqrt(1 + ei_abs2(u)); - rot1.s() = rot1.c() * u; - } - m.applyOnTheLeft(0,1,rot1); - j_right->makeJacobi(m,0,1); - *j_left = rot1 * j_right->transpose(); -} - template<typename MatrixType, int QRPreconditioner> JacobiSVD<MatrixType, QRPreconditioner>& JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsigned int computationOptions) @@ -551,8 +555,8 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig /*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */ - if(!ei_qr_preconditioner_impl<MatrixType, QRPreconditioner, PreconditionIfMoreColsThanRows>::run(*this, matrix) - && !ei_qr_preconditioner_impl<MatrixType, QRPreconditioner, PreconditionIfMoreRowsThanCols>::run(*this, matrix)) + if(!internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows>::run(*this, matrix) + && !internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols>::run(*this, matrix)) { m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize); if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows); @@ -577,15 +581,15 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig // if this 2x2 sub-matrix is not diagonal already... // notice that this comparison will evaluate to false if any NaN is involved, ensuring that NaN's don't // keep us iterating forever. - if(std::max(ei_abs(m_workMatrix.coeff(p,q)),ei_abs(m_workMatrix.coeff(q,p))) - > std::max(ei_abs(m_workMatrix.coeff(p,p)),ei_abs(m_workMatrix.coeff(q,q)))*precision) + if(std::max(internal::abs(m_workMatrix.coeff(p,q)),internal::abs(m_workMatrix.coeff(q,p))) + > std::max(internal::abs(m_workMatrix.coeff(p,p)),internal::abs(m_workMatrix.coeff(q,q)))*precision) { finished = false; // perform SVD decomposition of 2x2 sub-matrix corresponding to indices p,q to make it diagonal - ei_svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner>::run(m_workMatrix, *this, p, q); + internal::svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner>::run(m_workMatrix, *this, p, q); JacobiRotation<RealScalar> j_left, j_right; - ei_real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right); + internal::real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right); // accumulate resulting Jacobi rotations m_workMatrix.applyOnTheLeft(p,q,j_left); @@ -602,7 +606,7 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig for(Index i = 0; i < m_diagSize; ++i) { - RealScalar a = ei_abs(m_workMatrix.coeff(i,i)); + RealScalar a = internal::abs(m_workMatrix.coeff(i,i)); m_singularValues.coeffRef(i) = a; if(computeU() && (a!=RealScalar(0))) m_matrixU.col(i) *= m_workMatrix.coeff(i,i)/a; } @@ -632,16 +636,17 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig return *this; } +namespace internal { template<typename _MatrixType, int QRPreconditioner, typename Rhs> -struct ei_solve_retval<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs> - : ei_solve_retval_base<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs> +struct solve_retval<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs> + : solve_retval_base<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs> { typedef JacobiSVD<_MatrixType, QRPreconditioner> JacobiSVDType; EIGEN_MAKE_SOLVE_HELPERS(JacobiSVDType,Rhs) template<typename Dest> void evalTo(Dest& dst) const { - ei_assert(rhs().rows() == dec().rows()); + eigen_assert(rhs().rows() == dec().rows()); // A = U S V^* // So A^{-1} = V S^{-1} U^* @@ -659,6 +664,7 @@ struct ei_solve_retval<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs> * rhs(); } }; +} // end namespace internal template<typename Derived> JacobiSVD<typename MatrixBase<Derived>::PlainObject> |