FullPivLU.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2006-2009 Benoit Jacob <jacob.benoit.1@gmail.com>
00005 //
00006 // This Source Code Form is subject to the terms of the Mozilla
00007 // Public License v. 2.0. If a copy of the MPL was not distributed
00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00009 
00010 #ifndef EIGEN_LU_H
00011 #define EIGEN_LU_H
00012 
00013 namespace Eigen { 
00014 
00045 template<typename _MatrixType> class FullPivLU
00046 {
00047   public:
00048     typedef _MatrixType MatrixType;
00049     enum {
00050       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00051       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00052       Options = MatrixType::Options,
00053       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00054       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00055     };
00056     typedef typename MatrixType::Scalar Scalar;
00057     typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
00058     typedef typename internal::traits<MatrixType>::StorageKind StorageKind;
00059     typedef typename MatrixType::Index Index;
00060     typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType;
00061     typedef typename internal::plain_col_type<MatrixType, Index>::type IntColVectorType;
00062     typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationQType;
00063     typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationPType;
00064 
00071     FullPivLU();
00072 
00079     FullPivLU(Index rows, Index cols);
00080 
00086     FullPivLU(const MatrixType& matrix);
00087 
00095     FullPivLU& compute(const MatrixType& matrix);
00096 
00103     inline const MatrixType& matrixLU() const
00104     {
00105       eigen_assert(m_isInitialized && "LU is not initialized.");
00106       return m_lu;
00107     }
00108 
00116     inline Index nonzeroPivots() const
00117     {
00118       eigen_assert(m_isInitialized && "LU is not initialized.");
00119       return m_nonzero_pivots;
00120     }
00121 
00125     RealScalar maxPivot() const { return m_maxpivot; }
00126 
00131     inline const PermutationPType& permutationP() const
00132     {
00133       eigen_assert(m_isInitialized && "LU is not initialized.");
00134       return m_p;
00135     }
00136 
00141     inline const PermutationQType& permutationQ() const
00142     {
00143       eigen_assert(m_isInitialized && "LU is not initialized.");
00144       return m_q;
00145     }
00146 
00161     inline const internal::kernel_retval<FullPivLU> kernel() const
00162     {
00163       eigen_assert(m_isInitialized && "LU is not initialized.");
00164       return internal::kernel_retval<FullPivLU>(*this);
00165     }
00166 
00186     inline const internal::image_retval<FullPivLU>
00187       image(const MatrixType& originalMatrix) const
00188     {
00189       eigen_assert(m_isInitialized && "LU is not initialized.");
00190       return internal::image_retval<FullPivLU>(*this, originalMatrix);
00191     }
00192 
00212     template<typename Rhs>
00213     inline const internal::solve_retval<FullPivLU, Rhs>
00214     solve(const MatrixBase<Rhs>& b) const
00215     {
00216       eigen_assert(m_isInitialized && "LU is not initialized.");
00217       return internal::solve_retval<FullPivLU, Rhs>(*this, b.derived());
00218     }
00219 
00235     typename internal::traits<MatrixType>::Scalar determinant() const;
00236 
00254     FullPivLU& setThreshold(const RealScalar& threshold)
00255     {
00256       m_usePrescribedThreshold = true;
00257       m_prescribedThreshold = threshold;
00258       return *this;
00259     }
00260 
00269     FullPivLU& setThreshold(Default_t)
00270     {
00271       m_usePrescribedThreshold = false;
00272       return *this;
00273     }
00274 
00279     RealScalar threshold() const
00280     {
00281       eigen_assert(m_isInitialized || m_usePrescribedThreshold);
00282       return m_usePrescribedThreshold ? m_prescribedThreshold
00283       // this formula comes from experimenting (see "LU precision tuning" thread on the list)
00284       // and turns out to be identical to Higham's formula used already in LDLt.
00285                                       : NumTraits<Scalar>::epsilon() * m_lu.diagonalSize();
00286     }
00287 
00294     inline Index rank() const
00295     {
00296       using std::abs;
00297       eigen_assert(m_isInitialized && "LU is not initialized.");
00298       RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold();
00299       Index result = 0;
00300       for(Index i = 0; i < m_nonzero_pivots; ++i)
00301         result += (abs(m_lu.coeff(i,i)) > premultiplied_threshold);
00302       return result;
00303     }
00304 
00311     inline Index dimensionOfKernel() const
00312     {
00313       eigen_assert(m_isInitialized && "LU is not initialized.");
00314       return cols() - rank();
00315     }
00316 
00324     inline bool isInjective() const
00325     {
00326       eigen_assert(m_isInitialized && "LU is not initialized.");
00327       return rank() == cols();
00328     }
00329 
00337     inline bool isSurjective() const
00338     {
00339       eigen_assert(m_isInitialized && "LU is not initialized.");
00340       return rank() == rows();
00341     }
00342 
00349     inline bool isInvertible() const
00350     {
00351       eigen_assert(m_isInitialized && "LU is not initialized.");
00352       return isInjective() && (m_lu.rows() == m_lu.cols());
00353     }
00354 
00362     inline const internal::solve_retval<FullPivLU,typename MatrixType::IdentityReturnType> inverse() const
00363     {
00364       eigen_assert(m_isInitialized && "LU is not initialized.");
00365       eigen_assert(m_lu.rows() == m_lu.cols() && "You can't take the inverse of a non-square matrix!");
00366       return internal::solve_retval<FullPivLU,typename MatrixType::IdentityReturnType>
00367                (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols()));
00368     }
00369 
00370     MatrixType reconstructedMatrix() const;
00371 
00372     inline Index rows() const { return m_lu.rows(); }
00373     inline Index cols() const { return m_lu.cols(); }
00374 
00375   protected:
00376     MatrixType m_lu;
00377     PermutationPType m_p;
00378     PermutationQType m_q;
00379     IntColVectorType m_rowsTranspositions;
00380     IntRowVectorType m_colsTranspositions;
00381     Index m_det_pq, m_nonzero_pivots;
00382     RealScalar m_maxpivot, m_prescribedThreshold;
00383     bool m_isInitialized, m_usePrescribedThreshold;
00384 };
00385 
00386 template<typename MatrixType>
00387 FullPivLU<MatrixType>::FullPivLU()
00388   : m_isInitialized(false), m_usePrescribedThreshold(false)
00389 {
00390 }
00391 
00392 template<typename MatrixType>
00393 FullPivLU<MatrixType>::FullPivLU(Index rows, Index cols)
00394   : m_lu(rows, cols),
00395     m_p(rows),
00396     m_q(cols),
00397     m_rowsTranspositions(rows),
00398     m_colsTranspositions(cols),
00399     m_isInitialized(false),
00400     m_usePrescribedThreshold(false)
00401 {
00402 }
00403 
00404 template<typename MatrixType>
00405 FullPivLU<MatrixType>::FullPivLU(const MatrixType& matrix)
00406   : m_lu(matrix.rows(), matrix.cols()),
00407     m_p(matrix.rows()),
00408     m_q(matrix.cols()),
00409     m_rowsTranspositions(matrix.rows()),
00410     m_colsTranspositions(matrix.cols()),
00411     m_isInitialized(false),
00412     m_usePrescribedThreshold(false)
00413 {
00414   compute(matrix);
00415 }
00416 
00417 template<typename MatrixType>
00418 FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix)
00419 {
00420   // the permutations are stored as int indices, so just to be sure:
00421   eigen_assert(matrix.rows()<=NumTraits<int>::highest() && matrix.cols()<=NumTraits<int>::highest());
00422   
00423   m_isInitialized = true;
00424   m_lu = matrix;
00425 
00426   const Index size = matrix.diagonalSize();
00427   const Index rows = matrix.rows();
00428   const Index cols = matrix.cols();
00429 
00430   // will store the transpositions, before we accumulate them at the end.
00431   // can't accumulate on-the-fly because that will be done in reverse order for the rows.
00432   m_rowsTranspositions.resize(matrix.rows());
00433   m_colsTranspositions.resize(matrix.cols());
00434   Index number_of_transpositions = 0; // number of NONTRIVIAL transpositions, i.e. m_rowsTranspositions[i]!=i
00435 
00436   m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
00437   m_maxpivot = RealScalar(0);
00438 
00439   for(Index k = 0; k < size; ++k)
00440   {
00441     // First, we need to find the pivot.
00442 
00443     // biggest coefficient in the remaining bottom-right corner (starting at row k, col k)
00444     Index row_of_biggest_in_corner, col_of_biggest_in_corner;
00445     RealScalar biggest_in_corner;
00446     biggest_in_corner = m_lu.bottomRightCorner(rows-k, cols-k)
00447                         .cwiseAbs()
00448                         .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
00449     row_of_biggest_in_corner += k; // correct the values! since they were computed in the corner,
00450     col_of_biggest_in_corner += k; // need to add k to them.
00451 
00452     if(biggest_in_corner==RealScalar(0))
00453     {
00454       // before exiting, make sure to initialize the still uninitialized transpositions
00455       // in a sane state without destroying what we already have.
00456       m_nonzero_pivots = k;
00457       for(Index i = k; i < size; ++i)
00458       {
00459         m_rowsTranspositions.coeffRef(i) = i;
00460         m_colsTranspositions.coeffRef(i) = i;
00461       }
00462       break;
00463     }
00464 
00465     if(biggest_in_corner > m_maxpivot) m_maxpivot = biggest_in_corner;
00466 
00467     // Now that we've found the pivot, we need to apply the row/col swaps to
00468     // bring it to the location (k,k).
00469 
00470     m_rowsTranspositions.coeffRef(k) = row_of_biggest_in_corner;
00471     m_colsTranspositions.coeffRef(k) = col_of_biggest_in_corner;
00472     if(k != row_of_biggest_in_corner) {
00473       m_lu.row(k).swap(m_lu.row(row_of_biggest_in_corner));
00474       ++number_of_transpositions;
00475     }
00476     if(k != col_of_biggest_in_corner) {
00477       m_lu.col(k).swap(m_lu.col(col_of_biggest_in_corner));
00478       ++number_of_transpositions;
00479     }
00480 
00481     // Now that the pivot is at the right location, we update the remaining
00482     // bottom-right corner by Gaussian elimination.
00483 
00484     if(k<rows-1)
00485       m_lu.col(k).tail(rows-k-1) /= m_lu.coeff(k,k);
00486     if(k<size-1)
00487       m_lu.block(k+1,k+1,rows-k-1,cols-k-1).noalias() -= m_lu.col(k).tail(rows-k-1) * m_lu.row(k).tail(cols-k-1);
00488   }
00489 
00490   // the main loop is over, we still have to accumulate the transpositions to find the
00491   // permutations P and Q
00492 
00493   m_p.setIdentity(rows);
00494   for(Index k = size-1; k >= 0; --k)
00495     m_p.applyTranspositionOnTheRight(k, m_rowsTranspositions.coeff(k));
00496 
00497   m_q.setIdentity(cols);
00498   for(Index k = 0; k < size; ++k)
00499     m_q.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k));
00500 
00501   m_det_pq = (number_of_transpositions%2) ? -1 : 1;
00502   return *this;
00503 }
00504 
00505 template<typename MatrixType>
00506 typename internal::traits<MatrixType>::Scalar FullPivLU<MatrixType>::determinant() const
00507 {
00508   eigen_assert(m_isInitialized && "LU is not initialized.");
00509   eigen_assert(m_lu.rows() == m_lu.cols() && "You can't take the determinant of a non-square matrix!");
00510   return Scalar(m_det_pq) * Scalar(m_lu.diagonal().prod());
00511 }
00512 
00516 template<typename MatrixType>
00517 MatrixType FullPivLU<MatrixType>::reconstructedMatrix() const
00518 {
00519   eigen_assert(m_isInitialized && "LU is not initialized.");
00520   const Index smalldim = (std::min)(m_lu.rows(), m_lu.cols());
00521   // LU
00522   MatrixType res(m_lu.rows(),m_lu.cols());
00523   // FIXME the .toDenseMatrix() should not be needed...
00524   res = m_lu.leftCols(smalldim)
00525             .template triangularView<UnitLower>().toDenseMatrix()
00526       * m_lu.topRows(smalldim)
00527             .template triangularView<Upper>().toDenseMatrix();
00528 
00529   // P^{-1}(LU)
00530   res = m_p.inverse() * res;
00531 
00532   // (P^{-1}LU)Q^{-1}
00533   res = res * m_q.inverse();
00534 
00535   return res;
00536 }
00537 
00538 /********* Implementation of kernel() **************************************************/
00539 
00540 namespace internal {
00541 template<typename _MatrixType>
00542 struct kernel_retval<FullPivLU<_MatrixType> >
00543   : kernel_retval_base<FullPivLU<_MatrixType> >
00544 {
00545   EIGEN_MAKE_KERNEL_HELPERS(FullPivLU<_MatrixType>)
00546 
00547   enum { MaxSmallDimAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(
00548             MatrixType::MaxColsAtCompileTime,
00549             MatrixType::MaxRowsAtCompileTime)
00550   };
00551 
00552   template<typename Dest> void evalTo(Dest& dst) const
00553   {
00554     using std::abs;
00555     const Index cols = dec().matrixLU().cols(), dimker = cols - rank();
00556     if(dimker == 0)
00557     {
00558       // The Kernel is just {0}, so it doesn't have a basis properly speaking, but let's
00559       // avoid crashing/asserting as that depends on floating point calculations. Let's
00560       // just return a single column vector filled with zeros.
00561       dst.setZero();
00562       return;
00563     }
00564 
00565     /* Let us use the following lemma:
00566       *
00567       * Lemma: If the matrix A has the LU decomposition PAQ = LU,
00568       * then Ker A = Q(Ker U).
00569       *
00570       * Proof: trivial: just keep in mind that P, Q, L are invertible.
00571       */
00572 
00573     /* Thus, all we need to do is to compute Ker U, and then apply Q.
00574       *
00575       * U is upper triangular, with eigenvalues sorted so that any zeros appear at the end.
00576       * Thus, the diagonal of U ends with exactly
00577       * dimKer zero's. Let us use that to construct dimKer linearly
00578       * independent vectors in Ker U.
00579       */
00580 
00581     Matrix<Index, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank());
00582     RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold();
00583     Index p = 0;
00584     for(Index i = 0; i < dec().nonzeroPivots(); ++i)
00585       if(abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold)
00586         pivots.coeffRef(p++) = i;
00587     eigen_internal_assert(p == rank());
00588 
00589     // we construct a temporaty trapezoid matrix m, by taking the U matrix and
00590     // permuting the rows and cols to bring the nonnegligible pivots to the top of
00591     // the main diagonal. We need that to be able to apply our triangular solvers.
00592     // FIXME when we get triangularView-for-rectangular-matrices, this can be simplified
00593     Matrix<typename MatrixType::Scalar, Dynamic, Dynamic, MatrixType::Options,
00594            MaxSmallDimAtCompileTime, MatrixType::MaxColsAtCompileTime>
00595       m(dec().matrixLU().block(0, 0, rank(), cols));
00596     for(Index i = 0; i < rank(); ++i)
00597     {
00598       if(i) m.row(i).head(i).setZero();
00599       m.row(i).tail(cols-i) = dec().matrixLU().row(pivots.coeff(i)).tail(cols-i);
00600     }
00601     m.block(0, 0, rank(), rank());
00602     m.block(0, 0, rank(), rank()).template triangularView<StrictlyLower>().setZero();
00603     for(Index i = 0; i < rank(); ++i)
00604       m.col(i).swap(m.col(pivots.coeff(i)));
00605 
00606     // ok, we have our trapezoid matrix, we can apply the triangular solver.
00607     // notice that the math behind this suggests that we should apply this to the
00608     // negative of the RHS, but for performance we just put the negative sign elsewhere, see below.
00609     m.topLeftCorner(rank(), rank())
00610      .template triangularView<Upper>().solveInPlace(
00611         m.topRightCorner(rank(), dimker)
00612       );
00613 
00614     // now we must undo the column permutation that we had applied!
00615     for(Index i = rank()-1; i >= 0; --i)
00616       m.col(i).swap(m.col(pivots.coeff(i)));
00617 
00618     // see the negative sign in the next line, that's what we were talking about above.
00619     for(Index i = 0; i < rank(); ++i) dst.row(dec().permutationQ().indices().coeff(i)) = -m.row(i).tail(dimker);
00620     for(Index i = rank(); i < cols; ++i) dst.row(dec().permutationQ().indices().coeff(i)).setZero();
00621     for(Index k = 0; k < dimker; ++k) dst.coeffRef(dec().permutationQ().indices().coeff(rank()+k), k) = Scalar(1);
00622   }
00623 };
00624 
00625 /***** Implementation of image() *****************************************************/
00626 
00627 template<typename _MatrixType>
00628 struct image_retval<FullPivLU<_MatrixType> >
00629   : image_retval_base<FullPivLU<_MatrixType> >
00630 {
00631   EIGEN_MAKE_IMAGE_HELPERS(FullPivLU<_MatrixType>)
00632 
00633   enum { MaxSmallDimAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(
00634             MatrixType::MaxColsAtCompileTime,
00635             MatrixType::MaxRowsAtCompileTime)
00636   };
00637 
00638   template<typename Dest> void evalTo(Dest& dst) const
00639   {
00640     using std::abs;
00641     if(rank() == 0)
00642     {
00643       // The Image is just {0}, so it doesn't have a basis properly speaking, but let's
00644       // avoid crashing/asserting as that depends on floating point calculations. Let's
00645       // just return a single column vector filled with zeros.
00646       dst.setZero();
00647       return;
00648     }
00649 
00650     Matrix<Index, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank());
00651     RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold();
00652     Index p = 0;
00653     for(Index i = 0; i < dec().nonzeroPivots(); ++i)
00654       if(abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold)
00655         pivots.coeffRef(p++) = i;
00656     eigen_internal_assert(p == rank());
00657 
00658     for(Index i = 0; i < rank(); ++i)
00659       dst.col(i) = originalMatrix().col(dec().permutationQ().indices().coeff(pivots.coeff(i)));
00660   }
00661 };
00662 
00663 /***** Implementation of solve() *****************************************************/
00664 
00665 template<typename _MatrixType, typename Rhs>
00666 struct solve_retval<FullPivLU<_MatrixType>, Rhs>
00667   : solve_retval_base<FullPivLU<_MatrixType>, Rhs>
00668 {
00669   EIGEN_MAKE_SOLVE_HELPERS(FullPivLU<_MatrixType>,Rhs)
00670 
00671   template<typename Dest> void evalTo(Dest& dst) const
00672   {
00673     /* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1}.
00674      * So we proceed as follows:
00675      * Step 1: compute c = P * rhs.
00676      * Step 2: replace c by the solution x to Lx = c. Exists because L is invertible.
00677      * Step 3: replace c by the solution x to Ux = c. May or may not exist.
00678      * Step 4: result = Q * c;
00679      */
00680 
00681     const Index rows = dec().rows(), cols = dec().cols(),
00682               nonzero_pivots = dec().nonzeroPivots();
00683     eigen_assert(rhs().rows() == rows);
00684     const Index smalldim = (std::min)(rows, cols);
00685 
00686     if(nonzero_pivots == 0)
00687     {
00688       dst.setZero();
00689       return;
00690     }
00691 
00692     typename Rhs::PlainObject c(rhs().rows(), rhs().cols());
00693 
00694     // Step 1
00695     c = dec().permutationP() * rhs();
00696 
00697     // Step 2
00698     dec().matrixLU()
00699         .topLeftCorner(smalldim,smalldim)
00700         .template triangularView<UnitLower>()
00701         .solveInPlace(c.topRows(smalldim));
00702     if(rows>cols)
00703     {
00704       c.bottomRows(rows-cols)
00705         -= dec().matrixLU().bottomRows(rows-cols)
00706          * c.topRows(cols);
00707     }
00708 
00709     // Step 3
00710     dec().matrixLU()
00711         .topLeftCorner(nonzero_pivots, nonzero_pivots)
00712         .template triangularView<Upper>()
00713         .solveInPlace(c.topRows(nonzero_pivots));
00714 
00715     // Step 4
00716     for(Index i = 0; i < nonzero_pivots; ++i)
00717       dst.row(dec().permutationQ().indices().coeff(i)) = c.row(i);
00718     for(Index i = nonzero_pivots; i < dec().matrixLU().cols(); ++i)
00719       dst.row(dec().permutationQ().indices().coeff(i)).setZero();
00720   }
00721 };
00722 
00723 } // end namespace internal
00724 
00725 /******* MatrixBase methods *****************************************************************/
00726 
00733 template<typename Derived>
00734 inline const FullPivLU<typename MatrixBase<Derived>::PlainObject>
00735 MatrixBase<Derived>::fullPivLu() const
00736 {
00737   return FullPivLU<PlainObject>(eval());
00738 }
00739 
00740 } // end namespace Eigen
00741 
00742 #endif // EIGEN_LU_H


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Thu Aug 27 2015 11:58:19