Eigen  3.2.2
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
FullPivLU.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2006-2009 Benoit Jacob <jacob.benoit.1@gmail.com>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #ifndef EIGEN_LU_H
11 #define EIGEN_LU_H
12 
13 namespace Eigen {
14 
46 template<typename _MatrixType> class FullPivLU
47 {
48  public:
49  typedef _MatrixType MatrixType;
50  enum {
51  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
52  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
53  Options = MatrixType::Options,
54  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
55  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
56  };
57  typedef typename MatrixType::Scalar Scalar;
58  typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
59  typedef typename internal::traits<MatrixType>::StorageKind StorageKind;
60  typedef typename MatrixType::Index Index;
61  typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType;
62  typedef typename internal::plain_col_type<MatrixType, Index>::type IntColVectorType;
63  typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationQType;
64  typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationPType;
65 
72  FullPivLU();
73 
80  FullPivLU(Index rows, Index cols);
81 
87  FullPivLU(const MatrixType& matrix);
88 
96  FullPivLU& compute(const MatrixType& matrix);
97 
104  inline const MatrixType& matrixLU() const
105  {
106  eigen_assert(m_isInitialized && "LU is not initialized.");
107  return m_lu;
108  }
109 
117  inline Index nonzeroPivots() const
118  {
119  eigen_assert(m_isInitialized && "LU is not initialized.");
120  return m_nonzero_pivots;
121  }
122 
126  RealScalar maxPivot() const { return m_maxpivot; }
127 
132  inline const PermutationPType& permutationP() const
133  {
134  eigen_assert(m_isInitialized && "LU is not initialized.");
135  return m_p;
136  }
137 
142  inline const PermutationQType& permutationQ() const
143  {
144  eigen_assert(m_isInitialized && "LU is not initialized.");
145  return m_q;
146  }
147 
162  inline const internal::kernel_retval<FullPivLU> kernel() const
163  {
164  eigen_assert(m_isInitialized && "LU is not initialized.");
165  return internal::kernel_retval<FullPivLU>(*this);
166  }
167 
187  inline const internal::image_retval<FullPivLU>
188  image(const MatrixType& originalMatrix) const
189  {
190  eigen_assert(m_isInitialized && "LU is not initialized.");
191  return internal::image_retval<FullPivLU>(*this, originalMatrix);
192  }
193 
213  template<typename Rhs>
214  inline const internal::solve_retval<FullPivLU, Rhs>
215  solve(const MatrixBase<Rhs>& b) const
216  {
217  eigen_assert(m_isInitialized && "LU is not initialized.");
218  return internal::solve_retval<FullPivLU, Rhs>(*this, b.derived());
219  }
220 
236  typename internal::traits<MatrixType>::Scalar determinant() const;
237 
255  FullPivLU& setThreshold(const RealScalar& threshold)
256  {
257  m_usePrescribedThreshold = true;
258  m_prescribedThreshold = threshold;
259  return *this;
260  }
261 
271  {
272  m_usePrescribedThreshold = false;
273  return *this;
274  }
275 
280  RealScalar threshold() const
281  {
282  eigen_assert(m_isInitialized || m_usePrescribedThreshold);
283  return m_usePrescribedThreshold ? m_prescribedThreshold
284  // this formula comes from experimenting (see "LU precision tuning" thread on the list)
285  // and turns out to be identical to Higham's formula used already in LDLt.
286  : NumTraits<Scalar>::epsilon() * m_lu.diagonalSize();
287  }
288 
295  inline Index rank() const
296  {
297  using std::abs;
298  eigen_assert(m_isInitialized && "LU is not initialized.");
299  RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold();
300  Index result = 0;
301  for(Index i = 0; i < m_nonzero_pivots; ++i)
302  result += (abs(m_lu.coeff(i,i)) > premultiplied_threshold);
303  return result;
304  }
305 
312  inline Index dimensionOfKernel() const
313  {
314  eigen_assert(m_isInitialized && "LU is not initialized.");
315  return cols() - rank();
316  }
317 
325  inline bool isInjective() const
326  {
327  eigen_assert(m_isInitialized && "LU is not initialized.");
328  return rank() == cols();
329  }
330 
338  inline bool isSurjective() const
339  {
340  eigen_assert(m_isInitialized && "LU is not initialized.");
341  return rank() == rows();
342  }
343 
350  inline bool isInvertible() const
351  {
352  eigen_assert(m_isInitialized && "LU is not initialized.");
353  return isInjective() && (m_lu.rows() == m_lu.cols());
354  }
355 
363  inline const internal::solve_retval<FullPivLU,typename MatrixType::IdentityReturnType> inverse() const
364  {
365  eigen_assert(m_isInitialized && "LU is not initialized.");
366  eigen_assert(m_lu.rows() == m_lu.cols() && "You can't take the inverse of a non-square matrix!");
367  return internal::solve_retval<FullPivLU,typename MatrixType::IdentityReturnType>
368  (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols()));
369  }
370 
371  MatrixType reconstructedMatrix() const;
372 
373  inline Index rows() const { return m_lu.rows(); }
374  inline Index cols() const { return m_lu.cols(); }
375 
376  protected:
377  MatrixType m_lu;
378  PermutationPType m_p;
379  PermutationQType m_q;
380  IntColVectorType m_rowsTranspositions;
381  IntRowVectorType m_colsTranspositions;
382  Index m_det_pq, m_nonzero_pivots;
383  RealScalar m_maxpivot, m_prescribedThreshold;
384  bool m_isInitialized, m_usePrescribedThreshold;
385 };
386 
387 template<typename MatrixType>
389  : m_isInitialized(false), m_usePrescribedThreshold(false)
390 {
391 }
392 
393 template<typename MatrixType>
394 FullPivLU<MatrixType>::FullPivLU(Index rows, Index cols)
395  : m_lu(rows, cols),
396  m_p(rows),
397  m_q(cols),
398  m_rowsTranspositions(rows),
399  m_colsTranspositions(cols),
400  m_isInitialized(false),
401  m_usePrescribedThreshold(false)
402 {
403 }
404 
405 template<typename MatrixType>
406 FullPivLU<MatrixType>::FullPivLU(const MatrixType& matrix)
407  : m_lu(matrix.rows(), matrix.cols()),
408  m_p(matrix.rows()),
409  m_q(matrix.cols()),
410  m_rowsTranspositions(matrix.rows()),
411  m_colsTranspositions(matrix.cols()),
412  m_isInitialized(false),
413  m_usePrescribedThreshold(false)
414 {
415  compute(matrix);
416 }
417 
418 template<typename MatrixType>
420 {
421  // the permutations are stored as int indices, so just to be sure:
422  eigen_assert(matrix.rows()<=NumTraits<int>::highest() && matrix.cols()<=NumTraits<int>::highest());
423 
424  m_isInitialized = true;
425  m_lu = matrix;
426 
427  const Index size = matrix.diagonalSize();
428  const Index rows = matrix.rows();
429  const Index cols = matrix.cols();
430 
431  // will store the transpositions, before we accumulate them at the end.
432  // can't accumulate on-the-fly because that will be done in reverse order for the rows.
433  m_rowsTranspositions.resize(matrix.rows());
434  m_colsTranspositions.resize(matrix.cols());
435  Index number_of_transpositions = 0; // number of NONTRIVIAL transpositions, i.e. m_rowsTranspositions[i]!=i
436 
437  m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
438  m_maxpivot = RealScalar(0);
439 
440  for(Index k = 0; k < size; ++k)
441  {
442  // First, we need to find the pivot.
443 
444  // biggest coefficient in the remaining bottom-right corner (starting at row k, col k)
445  Index row_of_biggest_in_corner, col_of_biggest_in_corner;
446  RealScalar biggest_in_corner;
447  biggest_in_corner = m_lu.bottomRightCorner(rows-k, cols-k)
448  .cwiseAbs()
449  .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
450  row_of_biggest_in_corner += k; // correct the values! since they were computed in the corner,
451  col_of_biggest_in_corner += k; // need to add k to them.
452 
453  if(biggest_in_corner==RealScalar(0))
454  {
455  // before exiting, make sure to initialize the still uninitialized transpositions
456  // in a sane state without destroying what we already have.
457  m_nonzero_pivots = k;
458  for(Index i = k; i < size; ++i)
459  {
460  m_rowsTranspositions.coeffRef(i) = i;
461  m_colsTranspositions.coeffRef(i) = i;
462  }
463  break;
464  }
465 
466  if(biggest_in_corner > m_maxpivot) m_maxpivot = biggest_in_corner;
467 
468  // Now that we've found the pivot, we need to apply the row/col swaps to
469  // bring it to the location (k,k).
470 
471  m_rowsTranspositions.coeffRef(k) = row_of_biggest_in_corner;
472  m_colsTranspositions.coeffRef(k) = col_of_biggest_in_corner;
473  if(k != row_of_biggest_in_corner) {
474  m_lu.row(k).swap(m_lu.row(row_of_biggest_in_corner));
475  ++number_of_transpositions;
476  }
477  if(k != col_of_biggest_in_corner) {
478  m_lu.col(k).swap(m_lu.col(col_of_biggest_in_corner));
479  ++number_of_transpositions;
480  }
481 
482  // Now that the pivot is at the right location, we update the remaining
483  // bottom-right corner by Gaussian elimination.
484 
485  if(k<rows-1)
486  m_lu.col(k).tail(rows-k-1) /= m_lu.coeff(k,k);
487  if(k<size-1)
488  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);
489  }
490 
491  // the main loop is over, we still have to accumulate the transpositions to find the
492  // permutations P and Q
493 
494  m_p.setIdentity(rows);
495  for(Index k = size-1; k >= 0; --k)
496  m_p.applyTranspositionOnTheRight(k, m_rowsTranspositions.coeff(k));
497 
498  m_q.setIdentity(cols);
499  for(Index k = 0; k < size; ++k)
500  m_q.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k));
501 
502  m_det_pq = (number_of_transpositions%2) ? -1 : 1;
503  return *this;
504 }
505 
506 template<typename MatrixType>
507 typename internal::traits<MatrixType>::Scalar FullPivLU<MatrixType>::determinant() const
508 {
509  eigen_assert(m_isInitialized && "LU is not initialized.");
510  eigen_assert(m_lu.rows() == m_lu.cols() && "You can't take the determinant of a non-square matrix!");
511  return Scalar(m_det_pq) * Scalar(m_lu.diagonal().prod());
512 }
513 
517 template<typename MatrixType>
519 {
520  eigen_assert(m_isInitialized && "LU is not initialized.");
521  const Index smalldim = (std::min)(m_lu.rows(), m_lu.cols());
522  // LU
523  MatrixType res(m_lu.rows(),m_lu.cols());
524  // FIXME the .toDenseMatrix() should not be needed...
525  res = m_lu.leftCols(smalldim)
526  .template triangularView<UnitLower>().toDenseMatrix()
527  * m_lu.topRows(smalldim)
528  .template triangularView<Upper>().toDenseMatrix();
529 
530  // P^{-1}(LU)
531  res = m_p.inverse() * res;
532 
533  // (P^{-1}LU)Q^{-1}
534  res = res * m_q.inverse();
535 
536  return res;
537 }
538 
539 /********* Implementation of kernel() **************************************************/
540 
541 namespace internal {
542 template<typename _MatrixType>
543 struct kernel_retval<FullPivLU<_MatrixType> >
544  : kernel_retval_base<FullPivLU<_MatrixType> >
545 {
546  EIGEN_MAKE_KERNEL_HELPERS(FullPivLU<_MatrixType>)
547 
548  enum { MaxSmallDimAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(
549  MatrixType::MaxColsAtCompileTime,
550  MatrixType::MaxRowsAtCompileTime)
551  };
552 
553  template<typename Dest> void evalTo(Dest& dst) const
554  {
555  using std::abs;
556  const Index cols = dec().matrixLU().cols(), dimker = cols - rank();
557  if(dimker == 0)
558  {
559  // The Kernel is just {0}, so it doesn't have a basis properly speaking, but let's
560  // avoid crashing/asserting as that depends on floating point calculations. Let's
561  // just return a single column vector filled with zeros.
562  dst.setZero();
563  return;
564  }
565 
566  /* Let us use the following lemma:
567  *
568  * Lemma: If the matrix A has the LU decomposition PAQ = LU,
569  * then Ker A = Q(Ker U).
570  *
571  * Proof: trivial: just keep in mind that P, Q, L are invertible.
572  */
573 
574  /* Thus, all we need to do is to compute Ker U, and then apply Q.
575  *
576  * U is upper triangular, with eigenvalues sorted so that any zeros appear at the end.
577  * Thus, the diagonal of U ends with exactly
578  * dimKer zero's. Let us use that to construct dimKer linearly
579  * independent vectors in Ker U.
580  */
581 
582  Matrix<Index, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank());
583  RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold();
584  Index p = 0;
585  for(Index i = 0; i < dec().nonzeroPivots(); ++i)
586  if(abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold)
587  pivots.coeffRef(p++) = i;
588  eigen_internal_assert(p == rank());
589 
590  // we construct a temporaty trapezoid matrix m, by taking the U matrix and
591  // permuting the rows and cols to bring the nonnegligible pivots to the top of
592  // the main diagonal. We need that to be able to apply our triangular solvers.
593  // FIXME when we get triangularView-for-rectangular-matrices, this can be simplified
594  Matrix<typename MatrixType::Scalar, Dynamic, Dynamic, MatrixType::Options,
595  MaxSmallDimAtCompileTime, MatrixType::MaxColsAtCompileTime>
596  m(dec().matrixLU().block(0, 0, rank(), cols));
597  for(Index i = 0; i < rank(); ++i)
598  {
599  if(i) m.row(i).head(i).setZero();
600  m.row(i).tail(cols-i) = dec().matrixLU().row(pivots.coeff(i)).tail(cols-i);
601  }
602  m.block(0, 0, rank(), rank());
603  m.block(0, 0, rank(), rank()).template triangularView<StrictlyLower>().setZero();
604  for(Index i = 0; i < rank(); ++i)
605  m.col(i).swap(m.col(pivots.coeff(i)));
606 
607  // ok, we have our trapezoid matrix, we can apply the triangular solver.
608  // notice that the math behind this suggests that we should apply this to the
609  // negative of the RHS, but for performance we just put the negative sign elsewhere, see below.
610  m.topLeftCorner(rank(), rank())
611  .template triangularView<Upper>().solveInPlace(
612  m.topRightCorner(rank(), dimker)
613  );
614 
615  // now we must undo the column permutation that we had applied!
616  for(Index i = rank()-1; i >= 0; --i)
617  m.col(i).swap(m.col(pivots.coeff(i)));
618 
619  // see the negative sign in the next line, that's what we were talking about above.
620  for(Index i = 0; i < rank(); ++i) dst.row(dec().permutationQ().indices().coeff(i)) = -m.row(i).tail(dimker);
621  for(Index i = rank(); i < cols; ++i) dst.row(dec().permutationQ().indices().coeff(i)).setZero();
622  for(Index k = 0; k < dimker; ++k) dst.coeffRef(dec().permutationQ().indices().coeff(rank()+k), k) = Scalar(1);
623  }
624 };
625 
626 /***** Implementation of image() *****************************************************/
627 
628 template<typename _MatrixType>
629 struct image_retval<FullPivLU<_MatrixType> >
630  : image_retval_base<FullPivLU<_MatrixType> >
631 {
632  EIGEN_MAKE_IMAGE_HELPERS(FullPivLU<_MatrixType>)
633 
634  enum { MaxSmallDimAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(
635  MatrixType::MaxColsAtCompileTime,
636  MatrixType::MaxRowsAtCompileTime)
637  };
638 
639  template<typename Dest> void evalTo(Dest& dst) const
640  {
641  using std::abs;
642  if(rank() == 0)
643  {
644  // The Image is just {0}, so it doesn't have a basis properly speaking, but let's
645  // avoid crashing/asserting as that depends on floating point calculations. Let's
646  // just return a single column vector filled with zeros.
647  dst.setZero();
648  return;
649  }
650 
651  Matrix<Index, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank());
652  RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold();
653  Index p = 0;
654  for(Index i = 0; i < dec().nonzeroPivots(); ++i)
655  if(abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold)
656  pivots.coeffRef(p++) = i;
657  eigen_internal_assert(p == rank());
658 
659  for(Index i = 0; i < rank(); ++i)
660  dst.col(i) = originalMatrix().col(dec().permutationQ().indices().coeff(pivots.coeff(i)));
661  }
662 };
663 
664 /***** Implementation of solve() *****************************************************/
665 
666 template<typename _MatrixType, typename Rhs>
667 struct solve_retval<FullPivLU<_MatrixType>, Rhs>
668  : solve_retval_base<FullPivLU<_MatrixType>, Rhs>
669 {
670  EIGEN_MAKE_SOLVE_HELPERS(FullPivLU<_MatrixType>,Rhs)
671 
672  template<typename Dest> void evalTo(Dest& dst) const
673  {
674  /* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1}.
675  * So we proceed as follows:
676  * Step 1: compute c = P * rhs.
677  * Step 2: replace c by the solution x to Lx = c. Exists because L is invertible.
678  * Step 3: replace c by the solution x to Ux = c. May or may not exist.
679  * Step 4: result = Q * c;
680  */
681 
682  const Index rows = dec().rows(), cols = dec().cols(),
683  nonzero_pivots = dec().nonzeroPivots();
684  eigen_assert(rhs().rows() == rows);
685  const Index smalldim = (std::min)(rows, cols);
686 
687  if(nonzero_pivots == 0)
688  {
689  dst.setZero();
690  return;
691  }
692 
693  typename Rhs::PlainObject c(rhs().rows(), rhs().cols());
694 
695  // Step 1
696  c = dec().permutationP() * rhs();
697 
698  // Step 2
699  dec().matrixLU()
700  .topLeftCorner(smalldim,smalldim)
701  .template triangularView<UnitLower>()
702  .solveInPlace(c.topRows(smalldim));
703  if(rows>cols)
704  {
705  c.bottomRows(rows-cols)
706  -= dec().matrixLU().bottomRows(rows-cols)
707  * c.topRows(cols);
708  }
709 
710  // Step 3
711  dec().matrixLU()
712  .topLeftCorner(nonzero_pivots, nonzero_pivots)
713  .template triangularView<Upper>()
714  .solveInPlace(c.topRows(nonzero_pivots));
715 
716  // Step 4
717  for(Index i = 0; i < nonzero_pivots; ++i)
718  dst.row(dec().permutationQ().indices().coeff(i)) = c.row(i);
719  for(Index i = nonzero_pivots; i < dec().matrixLU().cols(); ++i)
720  dst.row(dec().permutationQ().indices().coeff(i)).setZero();
721  }
722 };
723 
724 } // end namespace internal
725 
726 /******* MatrixBase methods *****************************************************************/
727 
734 template<typename Derived>
735 inline const FullPivLU<typename MatrixBase<Derived>::PlainObject>
737 {
738  return FullPivLU<PlainObject>(eval());
739 }
740 
741 } // end namespace Eigen
742 
743 #endif // EIGEN_LU_H
bool isInvertible() const
Definition: FullPivLU.h:350
RealScalar threshold() const
Definition: FullPivLU.h:280
const MatrixType & matrixLU() const
Definition: FullPivLU.h:104
const FullPivLU< PlainObject > fullPivLu() const
Definition: FullPivLU.h:736
internal::traits< MatrixType >::Scalar determinant() const
Definition: FullPivLU.h:507
MatrixType reconstructedMatrix() const
Definition: FullPivLU.h:518
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:88
const int Dynamic
Definition: Constants.h:21
bool isInjective() const
Definition: FullPivLU.h:325
FullPivLU()
Default Constructor.
Definition: FullPivLU.h:388
FullPivLU & compute(const MatrixType &matrix)
Definition: FullPivLU.h:419
const PermutationPType & permutationP() const
Definition: FullPivLU.h:132
Index rank() const
Definition: FullPivLU.h:295
const internal::image_retval< FullPivLU > image(const MatrixType &originalMatrix) const
Definition: FullPivLU.h:188
FullPivLU & setThreshold(Default_t)
Definition: FullPivLU.h:270
Index nonzeroPivots() const
Definition: FullPivLU.h:117
Index dimensionOfKernel() const
Definition: FullPivLU.h:312
LU decomposition of a matrix with complete pivoting, and related features.
Definition: ForwardDeclarations.h:216
FullPivLU & setThreshold(const RealScalar &threshold)
Definition: FullPivLU.h:255
const internal::kernel_retval< FullPivLU > kernel() const
Definition: FullPivLU.h:162
RealScalar maxPivot() const
Definition: FullPivLU.h:126
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
const internal::solve_retval< FullPivLU, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: FullPivLU.h:215
bool isSurjective() const
Definition: FullPivLU.h:338
const internal::solve_retval< FullPivLU, typename MatrixType::IdentityReturnType > inverse() const
Definition: FullPivLU.h:363
const PermutationQType & permutationQ() const
Definition: FullPivLU.h:142