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 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #ifndef EIGEN_LU_H
00026 #define EIGEN_LU_H
00027 
00028 namespace Eigen { 
00029 
00060 template<typename _MatrixType> class FullPivLU
00061 {
00062   public:
00063     typedef _MatrixType MatrixType;
00064     enum {
00065       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00066       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00067       Options = MatrixType::Options,
00068       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00069       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00070     };
00071     typedef typename MatrixType::Scalar Scalar;
00072     typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
00073     typedef typename internal::traits<MatrixType>::StorageKind StorageKind;
00074     typedef typename MatrixType::Index Index;
00075     typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType;
00076     typedef typename internal::plain_col_type<MatrixType, Index>::type IntColVectorType;
00077     typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationQType;
00078     typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationPType;
00079 
00086     FullPivLU();
00087 
00094     FullPivLU(Index rows, Index cols);
00095 
00101     FullPivLU(const MatrixType& matrix);
00102 
00110     FullPivLU& compute(const MatrixType& matrix);
00111 
00118     inline const MatrixType& matrixLU() const
00119     {
00120       eigen_assert(m_isInitialized && "LU is not initialized.");
00121       return m_lu;
00122     }
00123 
00131     inline Index nonzeroPivots() const
00132     {
00133       eigen_assert(m_isInitialized && "LU is not initialized.");
00134       return m_nonzero_pivots;
00135     }
00136 
00140     RealScalar maxPivot() const { return m_maxpivot; }
00141 
00146     inline const PermutationPType& permutationP() const
00147     {
00148       eigen_assert(m_isInitialized && "LU is not initialized.");
00149       return m_p;
00150     }
00151 
00156     inline const PermutationQType& permutationQ() const
00157     {
00158       eigen_assert(m_isInitialized && "LU is not initialized.");
00159       return m_q;
00160     }
00161 
00176     inline const internal::kernel_retval<FullPivLU> kernel() const
00177     {
00178       eigen_assert(m_isInitialized && "LU is not initialized.");
00179       return internal::kernel_retval<FullPivLU>(*this);
00180     }
00181 
00201     inline const internal::image_retval<FullPivLU>
00202       image(const MatrixType& originalMatrix) const
00203     {
00204       eigen_assert(m_isInitialized && "LU is not initialized.");
00205       return internal::image_retval<FullPivLU>(*this, originalMatrix);
00206     }
00207 
00227     template<typename Rhs>
00228     inline const internal::solve_retval<FullPivLU, Rhs>
00229     solve(const MatrixBase<Rhs>& b) const
00230     {
00231       eigen_assert(m_isInitialized && "LU is not initialized.");
00232       return internal::solve_retval<FullPivLU, Rhs>(*this, b.derived());
00233     }
00234 
00250     typename internal::traits<MatrixType>::Scalar determinant() const;
00251 
00269     FullPivLU& setThreshold(const RealScalar& threshold)
00270     {
00271       m_usePrescribedThreshold = true;
00272       m_prescribedThreshold = threshold;
00273       return *this;
00274     }
00275 
00284     FullPivLU& setThreshold(Default_t)
00285     {
00286       m_usePrescribedThreshold = false;
00287       return *this;
00288     }
00289 
00294     RealScalar threshold() const
00295     {
00296       eigen_assert(m_isInitialized || m_usePrescribedThreshold);
00297       return m_usePrescribedThreshold ? m_prescribedThreshold
00298       // this formula comes from experimenting (see "LU precision tuning" thread on the list)
00299       // and turns out to be identical to Higham's formula used already in LDLt.
00300                                       : NumTraits<Scalar>::epsilon() * m_lu.diagonalSize();
00301     }
00302 
00309     inline Index rank() const
00310     {
00311       eigen_assert(m_isInitialized && "LU is not initialized.");
00312       RealScalar premultiplied_threshold = internal::abs(m_maxpivot) * threshold();
00313       Index result = 0;
00314       for(Index i = 0; i < m_nonzero_pivots; ++i)
00315         result += (internal::abs(m_lu.coeff(i,i)) > premultiplied_threshold);
00316       return result;
00317     }
00318 
00325     inline Index dimensionOfKernel() const
00326     {
00327       eigen_assert(m_isInitialized && "LU is not initialized.");
00328       return cols() - rank();
00329     }
00330 
00338     inline bool isInjective() const
00339     {
00340       eigen_assert(m_isInitialized && "LU is not initialized.");
00341       return rank() == cols();
00342     }
00343 
00351     inline bool isSurjective() const
00352     {
00353       eigen_assert(m_isInitialized && "LU is not initialized.");
00354       return rank() == rows();
00355     }
00356 
00363     inline bool isInvertible() const
00364     {
00365       eigen_assert(m_isInitialized && "LU is not initialized.");
00366       return isInjective() && (m_lu.rows() == m_lu.cols());
00367     }
00368 
00376     inline const internal::solve_retval<FullPivLU,typename MatrixType::IdentityReturnType> inverse() const
00377     {
00378       eigen_assert(m_isInitialized && "LU is not initialized.");
00379       eigen_assert(m_lu.rows() == m_lu.cols() && "You can't take the inverse of a non-square matrix!");
00380       return internal::solve_retval<FullPivLU,typename MatrixType::IdentityReturnType>
00381                (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols()));
00382     }
00383 
00384     MatrixType reconstructedMatrix() const;
00385 
00386     inline Index rows() const { return m_lu.rows(); }
00387     inline Index cols() const { return m_lu.cols(); }
00388 
00389   protected:
00390     MatrixType m_lu;
00391     PermutationPType m_p;
00392     PermutationQType m_q;
00393     IntColVectorType m_rowsTranspositions;
00394     IntRowVectorType m_colsTranspositions;
00395     Index m_det_pq, m_nonzero_pivots;
00396     RealScalar m_maxpivot, m_prescribedThreshold;
00397     bool m_isInitialized, m_usePrescribedThreshold;
00398 };
00399 
00400 template<typename MatrixType>
00401 FullPivLU<MatrixType>::FullPivLU()
00402   : m_isInitialized(false), m_usePrescribedThreshold(false)
00403 {
00404 }
00405 
00406 template<typename MatrixType>
00407 FullPivLU<MatrixType>::FullPivLU(Index rows, Index cols)
00408   : m_lu(rows, cols),
00409     m_p(rows),
00410     m_q(cols),
00411     m_rowsTranspositions(rows),
00412     m_colsTranspositions(cols),
00413     m_isInitialized(false),
00414     m_usePrescribedThreshold(false)
00415 {
00416 }
00417 
00418 template<typename MatrixType>
00419 FullPivLU<MatrixType>::FullPivLU(const MatrixType& matrix)
00420   : m_lu(matrix.rows(), matrix.cols()),
00421     m_p(matrix.rows()),
00422     m_q(matrix.cols()),
00423     m_rowsTranspositions(matrix.rows()),
00424     m_colsTranspositions(matrix.cols()),
00425     m_isInitialized(false),
00426     m_usePrescribedThreshold(false)
00427 {
00428   compute(matrix);
00429 }
00430 
00431 template<typename MatrixType>
00432 FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix)
00433 {
00434   m_isInitialized = true;
00435   m_lu = matrix;
00436 
00437   const Index size = matrix.diagonalSize();
00438   const Index rows = matrix.rows();
00439   const Index cols = matrix.cols();
00440 
00441   // will store the transpositions, before we accumulate them at the end.
00442   // can't accumulate on-the-fly because that will be done in reverse order for the rows.
00443   m_rowsTranspositions.resize(matrix.rows());
00444   m_colsTranspositions.resize(matrix.cols());
00445   Index number_of_transpositions = 0; // number of NONTRIVIAL transpositions, i.e. m_rowsTranspositions[i]!=i
00446 
00447   m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
00448   m_maxpivot = RealScalar(0);
00449 
00450   for(Index k = 0; k < size; ++k)
00451   {
00452     // First, we need to find the pivot.
00453 
00454     // biggest coefficient in the remaining bottom-right corner (starting at row k, col k)
00455     Index row_of_biggest_in_corner, col_of_biggest_in_corner;
00456     RealScalar biggest_in_corner;
00457     biggest_in_corner = m_lu.bottomRightCorner(rows-k, cols-k)
00458                         .cwiseAbs()
00459                         .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
00460     row_of_biggest_in_corner += k; // correct the values! since they were computed in the corner,
00461     col_of_biggest_in_corner += k; // need to add k to them.
00462 
00463     if(biggest_in_corner==RealScalar(0))
00464     {
00465       // before exiting, make sure to initialize the still uninitialized transpositions
00466       // in a sane state without destroying what we already have.
00467       m_nonzero_pivots = k;
00468       for(Index i = k; i < size; ++i)
00469       {
00470         m_rowsTranspositions.coeffRef(i) = i;
00471         m_colsTranspositions.coeffRef(i) = i;
00472       }
00473       break;
00474     }
00475 
00476     if(biggest_in_corner > m_maxpivot) m_maxpivot = biggest_in_corner;
00477 
00478     // Now that we've found the pivot, we need to apply the row/col swaps to
00479     // bring it to the location (k,k).
00480 
00481     m_rowsTranspositions.coeffRef(k) = row_of_biggest_in_corner;
00482     m_colsTranspositions.coeffRef(k) = col_of_biggest_in_corner;
00483     if(k != row_of_biggest_in_corner) {
00484       m_lu.row(k).swap(m_lu.row(row_of_biggest_in_corner));
00485       ++number_of_transpositions;
00486     }
00487     if(k != col_of_biggest_in_corner) {
00488       m_lu.col(k).swap(m_lu.col(col_of_biggest_in_corner));
00489       ++number_of_transpositions;
00490     }
00491 
00492     // Now that the pivot is at the right location, we update the remaining
00493     // bottom-right corner by Gaussian elimination.
00494 
00495     if(k<rows-1)
00496       m_lu.col(k).tail(rows-k-1) /= m_lu.coeff(k,k);
00497     if(k<size-1)
00498       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);
00499   }
00500 
00501   // the main loop is over, we still have to accumulate the transpositions to find the
00502   // permutations P and Q
00503 
00504   m_p.setIdentity(rows);
00505   for(Index k = size-1; k >= 0; --k)
00506     m_p.applyTranspositionOnTheRight(k, m_rowsTranspositions.coeff(k));
00507 
00508   m_q.setIdentity(cols);
00509   for(Index k = 0; k < size; ++k)
00510     m_q.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k));
00511 
00512   m_det_pq = (number_of_transpositions%2) ? -1 : 1;
00513   return *this;
00514 }
00515 
00516 template<typename MatrixType>
00517 typename internal::traits<MatrixType>::Scalar FullPivLU<MatrixType>::determinant() const
00518 {
00519   eigen_assert(m_isInitialized && "LU is not initialized.");
00520   eigen_assert(m_lu.rows() == m_lu.cols() && "You can't take the determinant of a non-square matrix!");
00521   return Scalar(m_det_pq) * Scalar(m_lu.diagonal().prod());
00522 }
00523 
00527 template<typename MatrixType>
00528 MatrixType FullPivLU<MatrixType>::reconstructedMatrix() const
00529 {
00530   eigen_assert(m_isInitialized && "LU is not initialized.");
00531   const Index smalldim = (std::min)(m_lu.rows(), m_lu.cols());
00532   // LU
00533   MatrixType res(m_lu.rows(),m_lu.cols());
00534   // FIXME the .toDenseMatrix() should not be needed...
00535   res = m_lu.leftCols(smalldim)
00536             .template triangularView<UnitLower>().toDenseMatrix()
00537       * m_lu.topRows(smalldim)
00538             .template triangularView<Upper>().toDenseMatrix();
00539 
00540   // P^{-1}(LU)
00541   res = m_p.inverse() * res;
00542 
00543   // (P^{-1}LU)Q^{-1}
00544   res = res * m_q.inverse();
00545 
00546   return res;
00547 }
00548 
00549 /********* Implementation of kernel() **************************************************/
00550 
00551 namespace internal {
00552 template<typename _MatrixType>
00553 struct kernel_retval<FullPivLU<_MatrixType> >
00554   : kernel_retval_base<FullPivLU<_MatrixType> >
00555 {
00556   EIGEN_MAKE_KERNEL_HELPERS(FullPivLU<_MatrixType>)
00557 
00558   enum { MaxSmallDimAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(
00559             MatrixType::MaxColsAtCompileTime,
00560             MatrixType::MaxRowsAtCompileTime)
00561   };
00562 
00563   template<typename Dest> void evalTo(Dest& dst) const
00564   {
00565     const Index cols = dec().matrixLU().cols(), dimker = cols - rank();
00566     if(dimker == 0)
00567     {
00568       // The Kernel is just {0}, so it doesn't have a basis properly speaking, but let's
00569       // avoid crashing/asserting as that depends on floating point calculations. Let's
00570       // just return a single column vector filled with zeros.
00571       dst.setZero();
00572       return;
00573     }
00574 
00575     /* Let us use the following lemma:
00576       *
00577       * Lemma: If the matrix A has the LU decomposition PAQ = LU,
00578       * then Ker A = Q(Ker U).
00579       *
00580       * Proof: trivial: just keep in mind that P, Q, L are invertible.
00581       */
00582 
00583     /* Thus, all we need to do is to compute Ker U, and then apply Q.
00584       *
00585       * U is upper triangular, with eigenvalues sorted so that any zeros appear at the end.
00586       * Thus, the diagonal of U ends with exactly
00587       * dimKer zero's. Let us use that to construct dimKer linearly
00588       * independent vectors in Ker U.
00589       */
00590 
00591     Matrix<Index, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank());
00592     RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold();
00593     Index p = 0;
00594     for(Index i = 0; i < dec().nonzeroPivots(); ++i)
00595       if(abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold)
00596         pivots.coeffRef(p++) = i;
00597     eigen_internal_assert(p == rank());
00598 
00599     // we construct a temporaty trapezoid matrix m, by taking the U matrix and
00600     // permuting the rows and cols to bring the nonnegligible pivots to the top of
00601     // the main diagonal. We need that to be able to apply our triangular solvers.
00602     // FIXME when we get triangularView-for-rectangular-matrices, this can be simplified
00603     Matrix<typename MatrixType::Scalar, Dynamic, Dynamic, MatrixType::Options,
00604            MaxSmallDimAtCompileTime, MatrixType::MaxColsAtCompileTime>
00605       m(dec().matrixLU().block(0, 0, rank(), cols));
00606     for(Index i = 0; i < rank(); ++i)
00607     {
00608       if(i) m.row(i).head(i).setZero();
00609       m.row(i).tail(cols-i) = dec().matrixLU().row(pivots.coeff(i)).tail(cols-i);
00610     }
00611     m.block(0, 0, rank(), rank());
00612     m.block(0, 0, rank(), rank()).template triangularView<StrictlyLower>().setZero();
00613     for(Index i = 0; i < rank(); ++i)
00614       m.col(i).swap(m.col(pivots.coeff(i)));
00615 
00616     // ok, we have our trapezoid matrix, we can apply the triangular solver.
00617     // notice that the math behind this suggests that we should apply this to the
00618     // negative of the RHS, but for performance we just put the negative sign elsewhere, see below.
00619     m.topLeftCorner(rank(), rank())
00620      .template triangularView<Upper>().solveInPlace(
00621         m.topRightCorner(rank(), dimker)
00622       );
00623 
00624     // now we must undo the column permutation that we had applied!
00625     for(Index i = rank()-1; i >= 0; --i)
00626       m.col(i).swap(m.col(pivots.coeff(i)));
00627 
00628     // see the negative sign in the next line, that's what we were talking about above.
00629     for(Index i = 0; i < rank(); ++i) dst.row(dec().permutationQ().indices().coeff(i)) = -m.row(i).tail(dimker);
00630     for(Index i = rank(); i < cols; ++i) dst.row(dec().permutationQ().indices().coeff(i)).setZero();
00631     for(Index k = 0; k < dimker; ++k) dst.coeffRef(dec().permutationQ().indices().coeff(rank()+k), k) = Scalar(1);
00632   }
00633 };
00634 
00635 /***** Implementation of image() *****************************************************/
00636 
00637 template<typename _MatrixType>
00638 struct image_retval<FullPivLU<_MatrixType> >
00639   : image_retval_base<FullPivLU<_MatrixType> >
00640 {
00641   EIGEN_MAKE_IMAGE_HELPERS(FullPivLU<_MatrixType>)
00642 
00643   enum { MaxSmallDimAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(
00644             MatrixType::MaxColsAtCompileTime,
00645             MatrixType::MaxRowsAtCompileTime)
00646   };
00647 
00648   template<typename Dest> void evalTo(Dest& dst) const
00649   {
00650     if(rank() == 0)
00651     {
00652       // The Image is just {0}, so it doesn't have a basis properly speaking, but let's
00653       // avoid crashing/asserting as that depends on floating point calculations. Let's
00654       // just return a single column vector filled with zeros.
00655       dst.setZero();
00656       return;
00657     }
00658 
00659     Matrix<Index, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank());
00660     RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold();
00661     Index p = 0;
00662     for(Index i = 0; i < dec().nonzeroPivots(); ++i)
00663       if(abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold)
00664         pivots.coeffRef(p++) = i;
00665     eigen_internal_assert(p == rank());
00666 
00667     for(Index i = 0; i < rank(); ++i)
00668       dst.col(i) = originalMatrix().col(dec().permutationQ().indices().coeff(pivots.coeff(i)));
00669   }
00670 };
00671 
00672 /***** Implementation of solve() *****************************************************/
00673 
00674 template<typename _MatrixType, typename Rhs>
00675 struct solve_retval<FullPivLU<_MatrixType>, Rhs>
00676   : solve_retval_base<FullPivLU<_MatrixType>, Rhs>
00677 {
00678   EIGEN_MAKE_SOLVE_HELPERS(FullPivLU<_MatrixType>,Rhs)
00679 
00680   template<typename Dest> void evalTo(Dest& dst) const
00681   {
00682     /* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1}.
00683      * So we proceed as follows:
00684      * Step 1: compute c = P * rhs.
00685      * Step 2: replace c by the solution x to Lx = c. Exists because L is invertible.
00686      * Step 3: replace c by the solution x to Ux = c. May or may not exist.
00687      * Step 4: result = Q * c;
00688      */
00689 
00690     const Index rows = dec().rows(), cols = dec().cols(),
00691               nonzero_pivots = dec().nonzeroPivots();
00692     eigen_assert(rhs().rows() == rows);
00693     const Index smalldim = (std::min)(rows, cols);
00694 
00695     if(nonzero_pivots == 0)
00696     {
00697       dst.setZero();
00698       return;
00699     }
00700 
00701     typename Rhs::PlainObject c(rhs().rows(), rhs().cols());
00702 
00703     // Step 1
00704     c = dec().permutationP() * rhs();
00705 
00706     // Step 2
00707     dec().matrixLU()
00708         .topLeftCorner(smalldim,smalldim)
00709         .template triangularView<UnitLower>()
00710         .solveInPlace(c.topRows(smalldim));
00711     if(rows>cols)
00712     {
00713       c.bottomRows(rows-cols)
00714         -= dec().matrixLU().bottomRows(rows-cols)
00715          * c.topRows(cols);
00716     }
00717 
00718     // Step 3
00719     dec().matrixLU()
00720         .topLeftCorner(nonzero_pivots, nonzero_pivots)
00721         .template triangularView<Upper>()
00722         .solveInPlace(c.topRows(nonzero_pivots));
00723 
00724     // Step 4
00725     for(Index i = 0; i < nonzero_pivots; ++i)
00726       dst.row(dec().permutationQ().indices().coeff(i)) = c.row(i);
00727     for(Index i = nonzero_pivots; i < dec().matrixLU().cols(); ++i)
00728       dst.row(dec().permutationQ().indices().coeff(i)).setZero();
00729   }
00730 };
00731 
00732 } // end namespace internal
00733 
00734 /******* MatrixBase methods *****************************************************************/
00735 
00742 template<typename Derived>
00743 inline const FullPivLU<typename MatrixBase<Derived>::PlainObject>
00744 MatrixBase<Derived>::fullPivLu() const
00745 {
00746   return FullPivLU<PlainObject>(eval());
00747 }
00748 
00749 } // end namespace Eigen
00750 
00751 #endif // EIGEN_LU_H