![]() |
Eigen
3.3.3
|
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2008-2011 Gael Guennebaud <gael.guennebaud@inria.fr> 00005 // Copyright (C) 2009 Keir Mierle <mierle@gmail.com> 00006 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> 00007 // Copyright (C) 2011 Timothy E. Holy <tim.holy@gmail.com > 00008 // 00009 // This Source Code Form is subject to the terms of the Mozilla 00010 // Public License v. 2.0. If a copy of the MPL was not distributed 00011 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00012 00013 #ifndef EIGEN_LDLT_H 00014 #define EIGEN_LDLT_H 00015 00016 namespace Eigen { 00017 00018 namespace internal { 00019 template<typename MatrixType, int UpLo> struct LDLT_Traits; 00020 00021 // PositiveSemiDef means positive semi-definite and non-zero; same for NegativeSemiDef 00022 enum SignMatrix { PositiveSemiDef, NegativeSemiDef, ZeroSign, Indefinite }; 00023 } 00024 00050 template<typename _MatrixType, int _UpLo> class LDLT 00051 { 00052 public: 00053 typedef _MatrixType MatrixType; 00054 enum { 00055 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 00056 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 00057 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 00058 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, 00059 UpLo = _UpLo 00060 }; 00061 typedef typename MatrixType::Scalar Scalar; 00062 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; 00063 typedef Eigen::Index Index; 00064 typedef typename MatrixType::StorageIndex StorageIndex; 00065 typedef Matrix<Scalar, RowsAtCompileTime, 1, 0, MaxRowsAtCompileTime, 1> TmpMatrixType; 00066 00067 typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType; 00068 typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType; 00069 00070 typedef internal::LDLT_Traits<MatrixType,UpLo> Traits; 00071 00077 LDLT() 00078 : m_matrix(), 00079 m_transpositions(), 00080 m_sign(internal::ZeroSign), 00081 m_isInitialized(false) 00082 {} 00083 00090 explicit LDLT(Index size) 00091 : m_matrix(size, size), 00092 m_transpositions(size), 00093 m_temporary(size), 00094 m_sign(internal::ZeroSign), 00095 m_isInitialized(false) 00096 {} 00097 00104 template<typename InputType> 00105 explicit LDLT(const EigenBase<InputType>& matrix) 00106 : m_matrix(matrix.rows(), matrix.cols()), 00107 m_transpositions(matrix.rows()), 00108 m_temporary(matrix.rows()), 00109 m_sign(internal::ZeroSign), 00110 m_isInitialized(false) 00111 { 00112 compute(matrix.derived()); 00113 } 00114 00121 template<typename InputType> 00122 explicit LDLT(EigenBase<InputType>& matrix) 00123 : m_matrix(matrix.derived()), 00124 m_transpositions(matrix.rows()), 00125 m_temporary(matrix.rows()), 00126 m_sign(internal::ZeroSign), 00127 m_isInitialized(false) 00128 { 00129 compute(matrix.derived()); 00130 } 00131 00135 void setZero() 00136 { 00137 m_isInitialized = false; 00138 } 00139 00141 inline typename Traits::MatrixU matrixU() const 00142 { 00143 eigen_assert(m_isInitialized && "LDLT is not initialized."); 00144 return Traits::getU(m_matrix); 00145 } 00146 00148 inline typename Traits::MatrixL matrixL() const 00149 { 00150 eigen_assert(m_isInitialized && "LDLT is not initialized."); 00151 return Traits::getL(m_matrix); 00152 } 00153 00156 inline const TranspositionType& transpositionsP() const 00157 { 00158 eigen_assert(m_isInitialized && "LDLT is not initialized."); 00159 return m_transpositions; 00160 } 00161 00163 inline Diagonal<const MatrixType> vectorD() const 00164 { 00165 eigen_assert(m_isInitialized && "LDLT is not initialized."); 00166 return m_matrix.diagonal(); 00167 } 00168 00170 inline bool isPositive() const 00171 { 00172 eigen_assert(m_isInitialized && "LDLT is not initialized."); 00173 return m_sign == internal::PositiveSemiDef || m_sign == internal::ZeroSign; 00174 } 00175 00177 inline bool isNegative(void) const 00178 { 00179 eigen_assert(m_isInitialized && "LDLT is not initialized."); 00180 return m_sign == internal::NegativeSemiDef || m_sign == internal::ZeroSign; 00181 } 00182 00198 template<typename Rhs> 00199 inline const Solve<LDLT, Rhs> 00200 solve(const MatrixBase<Rhs>& b) const 00201 { 00202 eigen_assert(m_isInitialized && "LDLT is not initialized."); 00203 eigen_assert(m_matrix.rows()==b.rows() 00204 && "LDLT::solve(): invalid number of rows of the right hand side matrix b"); 00205 return Solve<LDLT, Rhs>(*this, b.derived()); 00206 } 00207 00208 template<typename Derived> 00209 bool solveInPlace(MatrixBase<Derived> &bAndX) const; 00210 00211 template<typename InputType> 00212 LDLT& compute(const EigenBase<InputType>& matrix); 00213 00217 RealScalar rcond() const 00218 { 00219 eigen_assert(m_isInitialized && "LDLT is not initialized."); 00220 return internal::rcond_estimate_helper(m_l1_norm, *this); 00221 } 00222 00223 template <typename Derived> 00224 LDLT& rankUpdate(const MatrixBase<Derived>& w, const RealScalar& alpha=1); 00225 00230 inline const MatrixType& matrixLDLT() const 00231 { 00232 eigen_assert(m_isInitialized && "LDLT is not initialized."); 00233 return m_matrix; 00234 } 00235 00236 MatrixType reconstructedMatrix() const; 00237 00243 const LDLT& adjoint() const { return *this; }; 00244 00245 inline Index rows() const { return m_matrix.rows(); } 00246 inline Index cols() const { return m_matrix.cols(); } 00247 00253 ComputationInfo info() const 00254 { 00255 eigen_assert(m_isInitialized && "LDLT is not initialized."); 00256 return m_info; 00257 } 00258 00259 #ifndef EIGEN_PARSED_BY_DOXYGEN 00260 template<typename RhsType, typename DstType> 00261 EIGEN_DEVICE_FUNC 00262 void _solve_impl(const RhsType &rhs, DstType &dst) const; 00263 #endif 00264 00265 protected: 00266 00267 static void check_template_parameters() 00268 { 00269 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); 00270 } 00271 00278 MatrixType m_matrix; 00279 RealScalar m_l1_norm; 00280 TranspositionType m_transpositions; 00281 TmpMatrixType m_temporary; 00282 internal::SignMatrix m_sign; 00283 bool m_isInitialized; 00284 ComputationInfo m_info; 00285 }; 00286 00287 namespace internal { 00288 00289 template<int UpLo> struct ldlt_inplace; 00290 00291 template<> struct ldlt_inplace<Lower> 00292 { 00293 template<typename MatrixType, typename TranspositionType, typename Workspace> 00294 static bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, SignMatrix& sign) 00295 { 00296 using std::abs; 00297 typedef typename MatrixType::Scalar Scalar; 00298 typedef typename MatrixType::RealScalar RealScalar; 00299 typedef typename TranspositionType::StorageIndex IndexType; 00300 eigen_assert(mat.rows()==mat.cols()); 00301 const Index size = mat.rows(); 00302 bool found_zero_pivot = false; 00303 bool ret = true; 00304 00305 if (size <= 1) 00306 { 00307 transpositions.setIdentity(); 00308 if (numext::real(mat.coeff(0,0)) > static_cast<RealScalar>(0) ) sign = PositiveSemiDef; 00309 else if (numext::real(mat.coeff(0,0)) < static_cast<RealScalar>(0)) sign = NegativeSemiDef; 00310 else sign = ZeroSign; 00311 return true; 00312 } 00313 00314 for (Index k = 0; k < size; ++k) 00315 { 00316 // Find largest diagonal element 00317 Index index_of_biggest_in_corner; 00318 mat.diagonal().tail(size-k).cwiseAbs().maxCoeff(&index_of_biggest_in_corner); 00319 index_of_biggest_in_corner += k; 00320 00321 transpositions.coeffRef(k) = IndexType(index_of_biggest_in_corner); 00322 if(k != index_of_biggest_in_corner) 00323 { 00324 // apply the transposition while taking care to consider only 00325 // the lower triangular part 00326 Index s = size-index_of_biggest_in_corner-1; // trailing size after the biggest element 00327 mat.row(k).head(k).swap(mat.row(index_of_biggest_in_corner).head(k)); 00328 mat.col(k).tail(s).swap(mat.col(index_of_biggest_in_corner).tail(s)); 00329 std::swap(mat.coeffRef(k,k),mat.coeffRef(index_of_biggest_in_corner,index_of_biggest_in_corner)); 00330 for(Index i=k+1;i<index_of_biggest_in_corner;++i) 00331 { 00332 Scalar tmp = mat.coeffRef(i,k); 00333 mat.coeffRef(i,k) = numext::conj(mat.coeffRef(index_of_biggest_in_corner,i)); 00334 mat.coeffRef(index_of_biggest_in_corner,i) = numext::conj(tmp); 00335 } 00336 if(NumTraits<Scalar>::IsComplex) 00337 mat.coeffRef(index_of_biggest_in_corner,k) = numext::conj(mat.coeff(index_of_biggest_in_corner,k)); 00338 } 00339 00340 // partition the matrix: 00341 // A00 | - | - 00342 // lu = A10 | A11 | - 00343 // A20 | A21 | A22 00344 Index rs = size - k - 1; 00345 Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1); 00346 Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k); 00347 Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k); 00348 00349 if(k>0) 00350 { 00351 temp.head(k) = mat.diagonal().real().head(k).asDiagonal() * A10.adjoint(); 00352 mat.coeffRef(k,k) -= (A10 * temp.head(k)).value(); 00353 if(rs>0) 00354 A21.noalias() -= A20 * temp.head(k); 00355 } 00356 00357 // In some previous versions of Eigen (e.g., 3.2.1), the scaling was omitted if the pivot 00358 // was smaller than the cutoff value. However, since LDLT is not rank-revealing 00359 // we should only make sure that we do not introduce INF or NaN values. 00360 // Remark that LAPACK also uses 0 as the cutoff value. 00361 RealScalar realAkk = numext::real(mat.coeffRef(k,k)); 00362 bool pivot_is_valid = (abs(realAkk) > RealScalar(0)); 00363 00364 if(k==0 && !pivot_is_valid) 00365 { 00366 // The entire diagonal is zero, there is nothing more to do 00367 // except filling the transpositions, and checking whether the matrix is zero. 00368 sign = ZeroSign; 00369 for(Index j = 0; j<size; ++j) 00370 { 00371 transpositions.coeffRef(j) = IndexType(j); 00372 ret = ret && (mat.col(j).tail(size-j-1).array()==Scalar(0)).all(); 00373 } 00374 return ret; 00375 } 00376 00377 if((rs>0) && pivot_is_valid) 00378 A21 /= realAkk; 00379 00380 if(found_zero_pivot && pivot_is_valid) ret = false; // factorization failed 00381 else if(!pivot_is_valid) found_zero_pivot = true; 00382 00383 if (sign == PositiveSemiDef) { 00384 if (realAkk < static_cast<RealScalar>(0)) sign = Indefinite; 00385 } else if (sign == NegativeSemiDef) { 00386 if (realAkk > static_cast<RealScalar>(0)) sign = Indefinite; 00387 } else if (sign == ZeroSign) { 00388 if (realAkk > static_cast<RealScalar>(0)) sign = PositiveSemiDef; 00389 else if (realAkk < static_cast<RealScalar>(0)) sign = NegativeSemiDef; 00390 } 00391 } 00392 00393 return ret; 00394 } 00395 00396 // Reference for the algorithm: Davis and Hager, "Multiple Rank 00397 // Modifications of a Sparse Cholesky Factorization" (Algorithm 1) 00398 // Trivial rearrangements of their computations (Timothy E. Holy) 00399 // allow their algorithm to work for rank-1 updates even if the 00400 // original matrix is not of full rank. 00401 // Here only rank-1 updates are implemented, to reduce the 00402 // requirement for intermediate storage and improve accuracy 00403 template<typename MatrixType, typename WDerived> 00404 static bool updateInPlace(MatrixType& mat, MatrixBase<WDerived>& w, const typename MatrixType::RealScalar& sigma=1) 00405 { 00406 using numext::isfinite; 00407 typedef typename MatrixType::Scalar Scalar; 00408 typedef typename MatrixType::RealScalar RealScalar; 00409 00410 const Index size = mat.rows(); 00411 eigen_assert(mat.cols() == size && w.size()==size); 00412 00413 RealScalar alpha = 1; 00414 00415 // Apply the update 00416 for (Index j = 0; j < size; j++) 00417 { 00418 // Check for termination due to an original decomposition of low-rank 00419 if (!(isfinite)(alpha)) 00420 break; 00421 00422 // Update the diagonal terms 00423 RealScalar dj = numext::real(mat.coeff(j,j)); 00424 Scalar wj = w.coeff(j); 00425 RealScalar swj2 = sigma*numext::abs2(wj); 00426 RealScalar gamma = dj*alpha + swj2; 00427 00428 mat.coeffRef(j,j) += swj2/alpha; 00429 alpha += swj2/dj; 00430 00431 00432 // Update the terms of L 00433 Index rs = size-j-1; 00434 w.tail(rs) -= wj * mat.col(j).tail(rs); 00435 if(gamma != 0) 00436 mat.col(j).tail(rs) += (sigma*numext::conj(wj)/gamma)*w.tail(rs); 00437 } 00438 return true; 00439 } 00440 00441 template<typename MatrixType, typename TranspositionType, typename Workspace, typename WType> 00442 static bool update(MatrixType& mat, const TranspositionType& transpositions, Workspace& tmp, const WType& w, const typename MatrixType::RealScalar& sigma=1) 00443 { 00444 // Apply the permutation to the input w 00445 tmp = transpositions * w; 00446 00447 return ldlt_inplace<Lower>::updateInPlace(mat,tmp,sigma); 00448 } 00449 }; 00450 00451 template<> struct ldlt_inplace<Upper> 00452 { 00453 template<typename MatrixType, typename TranspositionType, typename Workspace> 00454 static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, SignMatrix& sign) 00455 { 00456 Transpose<MatrixType> matt(mat); 00457 return ldlt_inplace<Lower>::unblocked(matt, transpositions, temp, sign); 00458 } 00459 00460 template<typename MatrixType, typename TranspositionType, typename Workspace, typename WType> 00461 static EIGEN_STRONG_INLINE bool update(MatrixType& mat, TranspositionType& transpositions, Workspace& tmp, WType& w, const typename MatrixType::RealScalar& sigma=1) 00462 { 00463 Transpose<MatrixType> matt(mat); 00464 return ldlt_inplace<Lower>::update(matt, transpositions, tmp, w.conjugate(), sigma); 00465 } 00466 }; 00467 00468 template<typename MatrixType> struct LDLT_Traits<MatrixType,Lower> 00469 { 00470 typedef const TriangularView<const MatrixType, UnitLower> MatrixL; 00471 typedef const TriangularView<const typename MatrixType::AdjointReturnType, UnitUpper> MatrixU; 00472 static inline MatrixL getL(const MatrixType& m) { return MatrixL(m); } 00473 static inline MatrixU getU(const MatrixType& m) { return MatrixU(m.adjoint()); } 00474 }; 00475 00476 template<typename MatrixType> struct LDLT_Traits<MatrixType,Upper> 00477 { 00478 typedef const TriangularView<const typename MatrixType::AdjointReturnType, UnitLower> MatrixL; 00479 typedef const TriangularView<const MatrixType, UnitUpper> MatrixU; 00480 static inline MatrixL getL(const MatrixType& m) { return MatrixL(m.adjoint()); } 00481 static inline MatrixU getU(const MatrixType& m) { return MatrixU(m); } 00482 }; 00483 00484 } // end namespace internal 00485 00488 template<typename MatrixType, int _UpLo> 00489 template<typename InputType> 00490 LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::compute(const EigenBase<InputType>& a) 00491 { 00492 check_template_parameters(); 00493 00494 eigen_assert(a.rows()==a.cols()); 00495 const Index size = a.rows(); 00496 00497 m_matrix = a.derived(); 00498 00499 // Compute matrix L1 norm = max abs column sum. 00500 m_l1_norm = RealScalar(0); 00501 // TODO move this code to SelfAdjointView 00502 for (Index col = 0; col < size; ++col) { 00503 RealScalar abs_col_sum; 00504 if (_UpLo == Lower) 00505 abs_col_sum = m_matrix.col(col).tail(size - col).template lpNorm<1>() + m_matrix.row(col).head(col).template lpNorm<1>(); 00506 else 00507 abs_col_sum = m_matrix.col(col).head(col).template lpNorm<1>() + m_matrix.row(col).tail(size - col).template lpNorm<1>(); 00508 if (abs_col_sum > m_l1_norm) 00509 m_l1_norm = abs_col_sum; 00510 } 00511 00512 m_transpositions.resize(size); 00513 m_isInitialized = false; 00514 m_temporary.resize(size); 00515 m_sign = internal::ZeroSign; 00516 00517 m_info = internal::ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, m_sign) ? Success : NumericalIssue; 00518 00519 m_isInitialized = true; 00520 return *this; 00521 } 00522 00528 template<typename MatrixType, int _UpLo> 00529 template<typename Derived> 00530 LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::rankUpdate(const MatrixBase<Derived>& w, const typename LDLT<MatrixType,_UpLo>::RealScalar& sigma) 00531 { 00532 typedef typename TranspositionType::StorageIndex IndexType; 00533 const Index size = w.rows(); 00534 if (m_isInitialized) 00535 { 00536 eigen_assert(m_matrix.rows()==size); 00537 } 00538 else 00539 { 00540 m_matrix.resize(size,size); 00541 m_matrix.setZero(); 00542 m_transpositions.resize(size); 00543 for (Index i = 0; i < size; i++) 00544 m_transpositions.coeffRef(i) = IndexType(i); 00545 m_temporary.resize(size); 00546 m_sign = sigma>=0 ? internal::PositiveSemiDef : internal::NegativeSemiDef; 00547 m_isInitialized = true; 00548 } 00549 00550 internal::ldlt_inplace<UpLo>::update(m_matrix, m_transpositions, m_temporary, w, sigma); 00551 00552 return *this; 00553 } 00554 00555 #ifndef EIGEN_PARSED_BY_DOXYGEN 00556 template<typename _MatrixType, int _UpLo> 00557 template<typename RhsType, typename DstType> 00558 void LDLT<_MatrixType,_UpLo>::_solve_impl(const RhsType &rhs, DstType &dst) const 00559 { 00560 eigen_assert(rhs.rows() == rows()); 00561 // dst = P b 00562 dst = m_transpositions * rhs; 00563 00564 // dst = L^-1 (P b) 00565 matrixL().solveInPlace(dst); 00566 00567 // dst = D^-1 (L^-1 P b) 00568 // more precisely, use pseudo-inverse of D (see bug 241) 00569 using std::abs; 00570 const typename Diagonal<const MatrixType>::RealReturnType vecD(vectorD()); 00571 // In some previous versions, tolerance was set to the max of 1/highest and the maximal diagonal entry * epsilon 00572 // as motivated by LAPACK's xGELSS: 00573 // RealScalar tolerance = numext::maxi(vecD.array().abs().maxCoeff() * NumTraits<RealScalar>::epsilon(),RealScalar(1) / NumTraits<RealScalar>::highest()); 00574 // However, LDLT is not rank revealing, and so adjusting the tolerance wrt to the highest 00575 // diagonal element is not well justified and leads to numerical issues in some cases. 00576 // Moreover, Lapack's xSYTRS routines use 0 for the tolerance. 00577 RealScalar tolerance = RealScalar(1) / NumTraits<RealScalar>::highest(); 00578 00579 for (Index i = 0; i < vecD.size(); ++i) 00580 { 00581 if(abs(vecD(i)) > tolerance) 00582 dst.row(i) /= vecD(i); 00583 else 00584 dst.row(i).setZero(); 00585 } 00586 00587 // dst = L^-T (D^-1 L^-1 P b) 00588 matrixU().solveInPlace(dst); 00589 00590 // dst = P^-1 (L^-T D^-1 L^-1 P b) = A^-1 b 00591 dst = m_transpositions.transpose() * dst; 00592 } 00593 #endif 00594 00608 template<typename MatrixType,int _UpLo> 00609 template<typename Derived> 00610 bool LDLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const 00611 { 00612 eigen_assert(m_isInitialized && "LDLT is not initialized."); 00613 eigen_assert(m_matrix.rows() == bAndX.rows()); 00614 00615 bAndX = this->solve(bAndX); 00616 00617 return true; 00618 } 00619 00623 template<typename MatrixType, int _UpLo> 00624 MatrixType LDLT<MatrixType,_UpLo>::reconstructedMatrix() const 00625 { 00626 eigen_assert(m_isInitialized && "LDLT is not initialized."); 00627 const Index size = m_matrix.rows(); 00628 MatrixType res(size,size); 00629 00630 // P 00631 res.setIdentity(); 00632 res = transpositionsP() * res; 00633 // L^* P 00634 res = matrixU() * res; 00635 // D(L^*P) 00636 res = vectorD().real().asDiagonal() * res; 00637 // L(DL^*P) 00638 res = matrixL() * res; 00639 // P^T (LDL^*P) 00640 res = transpositionsP().transpose() * res; 00641 00642 return res; 00643 } 00644 00649 template<typename MatrixType, unsigned int UpLo> 00650 inline const LDLT<typename SelfAdjointView<MatrixType, UpLo>::PlainObject, UpLo> 00651 SelfAdjointView<MatrixType, UpLo>::ldlt() const 00652 { 00653 return LDLT<PlainObject,UpLo>(m_matrix); 00654 } 00655 00660 template<typename Derived> 00661 inline const LDLT<typename MatrixBase<Derived>::PlainObject> 00662 MatrixBase<Derived>::ldlt() const 00663 { 00664 return LDLT<PlainObject>(derived()); 00665 } 00666 00667 } // end namespace Eigen 00668 00669 #endif // EIGEN_LDLT_H