FullPivLU.h
Go to the documentation of this file.
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 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24 
25 #ifndef EIGEN_LU_H
26 #define EIGEN_LU_H
27 
28 namespace Eigen {
29 
60 template<typename _MatrixType> class FullPivLU
61 {
62  public:
63  typedef _MatrixType MatrixType;
64  enum {
69  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
70  };
71  typedef typename MatrixType::Scalar Scalar;
73  typedef typename internal::traits<MatrixType>::StorageKind StorageKind;
74  typedef typename MatrixType::Index Index;
75  typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType;
76  typedef typename internal::plain_col_type<MatrixType, Index>::type IntColVectorType;
79 
86  FullPivLU();
87 
95 
101  FullPivLU(const MatrixType& matrix);
102 
110  FullPivLU& compute(const MatrixType& matrix);
111 
118  inline const MatrixType& matrixLU() const
119  {
120  eigen_assert(m_isInitialized && "LU is not initialized.");
121  return m_lu;
122  }
123 
131  inline Index nonzeroPivots() const
132  {
133  eigen_assert(m_isInitialized && "LU is not initialized.");
134  return m_nonzero_pivots;
135  }
136 
140  RealScalar maxPivot() const { return m_maxpivot; }
141 
146  inline const PermutationPType& permutationP() const
147  {
148  eigen_assert(m_isInitialized && "LU is not initialized.");
149  return m_p;
150  }
151 
156  inline const PermutationQType& permutationQ() const
157  {
158  eigen_assert(m_isInitialized && "LU is not initialized.");
159  return m_q;
160  }
161 
176  inline const internal::kernel_retval<FullPivLU> kernel() const
177  {
178  eigen_assert(m_isInitialized && "LU is not initialized.");
179  return internal::kernel_retval<FullPivLU>(*this);
180  }
181 
201  inline const internal::image_retval<FullPivLU>
202  image(const MatrixType& originalMatrix) const
203  {
204  eigen_assert(m_isInitialized && "LU is not initialized.");
205  return internal::image_retval<FullPivLU>(*this, originalMatrix);
206  }
207 
227  template<typename Rhs>
228  inline const internal::solve_retval<FullPivLU, Rhs>
229  solve(const MatrixBase<Rhs>& b) const
230  {
231  eigen_assert(m_isInitialized && "LU is not initialized.");
232  return internal::solve_retval<FullPivLU, Rhs>(*this, b.derived());
233  }
234 
250  typename internal::traits<MatrixType>::Scalar determinant() const;
251 
270  {
273  return *this;
274  }
275 
285  {
286  m_usePrescribedThreshold = false;
287  return *this;
288  }
289 
295  {
298  // this formula comes from experimenting (see "LU precision tuning" thread on the list)
299  // and turns out to be identical to Higham's formula used already in LDLt.
300  : NumTraits<Scalar>::epsilon() * m_lu.diagonalSize();
301  }
302 
309  inline Index rank() const
310  {
311  eigen_assert(m_isInitialized && "LU is not initialized.");
312  RealScalar premultiplied_threshold = internal::abs(m_maxpivot) * threshold();
313  Index result = 0;
314  for(Index i = 0; i < m_nonzero_pivots; ++i)
315  result += (internal::abs(m_lu.coeff(i,i)) > premultiplied_threshold);
316  return result;
317  }
318 
325  inline Index dimensionOfKernel() const
326  {
327  eigen_assert(m_isInitialized && "LU is not initialized.");
328  return cols() - rank();
329  }
330 
338  inline bool isInjective() const
339  {
340  eigen_assert(m_isInitialized && "LU is not initialized.");
341  return rank() == cols();
342  }
343 
351  inline bool isSurjective() const
352  {
353  eigen_assert(m_isInitialized && "LU is not initialized.");
354  return rank() == rows();
355  }
356 
363  inline bool isInvertible() const
364  {
365  eigen_assert(m_isInitialized && "LU is not initialized.");
366  return isInjective() && (m_lu.rows() == m_lu.cols());
367  }
368 
376  inline const internal::solve_retval<FullPivLU,typename MatrixType::IdentityReturnType> inverse() const
377  {
378  eigen_assert(m_isInitialized && "LU is not initialized.");
379  eigen_assert(m_lu.rows() == m_lu.cols() && "You can't take the inverse of a non-square matrix!");
380  return internal::solve_retval<FullPivLU,typename MatrixType::IdentityReturnType>
381  (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols()));
382  }
383 
385 
386  inline Index rows() const { return m_lu.rows(); }
387  inline Index cols() const { return m_lu.cols(); }
388 
389  protected:
398 };
399 
400 template<typename MatrixType>
402  : m_isInitialized(false), m_usePrescribedThreshold(false)
403 {
404 }
405 
406 template<typename MatrixType>
408  : m_lu(rows, cols),
409  m_p(rows),
410  m_q(cols),
411  m_rowsTranspositions(rows),
412  m_colsTranspositions(cols),
413  m_isInitialized(false),
414  m_usePrescribedThreshold(false)
415 {
416 }
417 
418 template<typename MatrixType>
420  : m_lu(matrix.rows(), matrix.cols()),
421  m_p(matrix.rows()),
422  m_q(matrix.cols()),
423  m_rowsTranspositions(matrix.rows()),
424  m_colsTranspositions(matrix.cols()),
425  m_isInitialized(false),
426  m_usePrescribedThreshold(false)
427 {
428  compute(matrix);
429 }
430 
431 template<typename MatrixType>
433 {
434  m_isInitialized = true;
435  m_lu = matrix;
436 
437  const Index size = matrix.diagonalSize();
438  const Index rows = matrix.rows();
439  const Index cols = matrix.cols();
440 
441  // will store the transpositions, before we accumulate them at the end.
442  // can't accumulate on-the-fly because that will be done in reverse order for the rows.
443  m_rowsTranspositions.resize(matrix.rows());
444  m_colsTranspositions.resize(matrix.cols());
445  Index number_of_transpositions = 0; // number of NONTRIVIAL transpositions, i.e. m_rowsTranspositions[i]!=i
446 
447  m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
448  m_maxpivot = RealScalar(0);
449 
450  for(Index k = 0; k < size; ++k)
451  {
452  // First, we need to find the pivot.
453 
454  // biggest coefficient in the remaining bottom-right corner (starting at row k, col k)
455  Index row_of_biggest_in_corner, col_of_biggest_in_corner;
456  RealScalar biggest_in_corner;
457  biggest_in_corner = m_lu.bottomRightCorner(rows-k, cols-k)
458  .cwiseAbs()
459  .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
460  row_of_biggest_in_corner += k; // correct the values! since they were computed in the corner,
461  col_of_biggest_in_corner += k; // need to add k to them.
462 
463  if(biggest_in_corner==RealScalar(0))
464  {
465  // before exiting, make sure to initialize the still uninitialized transpositions
466  // in a sane state without destroying what we already have.
467  m_nonzero_pivots = k;
468  for(Index i = k; i < size; ++i)
469  {
470  m_rowsTranspositions.coeffRef(i) = i;
471  m_colsTranspositions.coeffRef(i) = i;
472  }
473  break;
474  }
475 
476  if(biggest_in_corner > m_maxpivot) m_maxpivot = biggest_in_corner;
477 
478  // Now that we've found the pivot, we need to apply the row/col swaps to
479  // bring it to the location (k,k).
480 
481  m_rowsTranspositions.coeffRef(k) = row_of_biggest_in_corner;
482  m_colsTranspositions.coeffRef(k) = col_of_biggest_in_corner;
483  if(k != row_of_biggest_in_corner) {
484  m_lu.row(k).swap(m_lu.row(row_of_biggest_in_corner));
485  ++number_of_transpositions;
486  }
487  if(k != col_of_biggest_in_corner) {
488  m_lu.col(k).swap(m_lu.col(col_of_biggest_in_corner));
489  ++number_of_transpositions;
490  }
491 
492  // Now that the pivot is at the right location, we update the remaining
493  // bottom-right corner by Gaussian elimination.
494 
495  if(k<rows-1)
496  m_lu.col(k).tail(rows-k-1) /= m_lu.coeff(k,k);
497  if(k<size-1)
498  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);
499  }
500 
501  // the main loop is over, we still have to accumulate the transpositions to find the
502  // permutations P and Q
503 
504  m_p.setIdentity(rows);
505  for(Index k = size-1; k >= 0; --k)
506  m_p.applyTranspositionOnTheRight(k, m_rowsTranspositions.coeff(k));
507 
508  m_q.setIdentity(cols);
509  for(Index k = 0; k < size; ++k)
510  m_q.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k));
511 
512  m_det_pq = (number_of_transpositions%2) ? -1 : 1;
513  return *this;
514 }
515 
516 template<typename MatrixType>
517 typename internal::traits<MatrixType>::Scalar FullPivLU<MatrixType>::determinant() const
518 {
519  eigen_assert(m_isInitialized && "LU is not initialized.");
520  eigen_assert(m_lu.rows() == m_lu.cols() && "You can't take the determinant of a non-square matrix!");
521  return Scalar(m_det_pq) * Scalar(m_lu.diagonal().prod());
522 }
523 
527 template<typename MatrixType>
529 {
530  eigen_assert(m_isInitialized && "LU is not initialized.");
531  const Index smalldim = (std::min)(m_lu.rows(), m_lu.cols());
532  // LU
533  MatrixType res(m_lu.rows(),m_lu.cols());
534  // FIXME the .toDenseMatrix() should not be needed...
535  res = m_lu.leftCols(smalldim)
536  .template triangularView<UnitLower>().toDenseMatrix()
537  * m_lu.topRows(smalldim)
538  .template triangularView<Upper>().toDenseMatrix();
539 
540  // P^{-1}(LU)
541  res = m_p.inverse() * res;
542 
543  // (P^{-1}LU)Q^{-1}
544  res = res * m_q.inverse();
545 
546  return res;
547 }
548 
549 /********* Implementation of kernel() **************************************************/
550 
551 namespace internal {
552 template<typename _MatrixType>
553 struct kernel_retval<FullPivLU<_MatrixType> >
554  : kernel_retval_base<FullPivLU<_MatrixType> >
555 {
557 
558  enum { MaxSmallDimAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(
559  MatrixType::MaxColsAtCompileTime,
560  MatrixType::MaxRowsAtCompileTime)
561  };
562 
563  template<typename Dest> void evalTo(Dest& dst) const
564  {
565  const Index cols = dec().matrixLU().cols(), dimker = cols - rank();
566  if(dimker == 0)
567  {
568  // The Kernel is just {0}, so it doesn't have a basis properly speaking, but let's
569  // avoid crashing/asserting as that depends on floating point calculations. Let's
570  // just return a single column vector filled with zeros.
571  dst.setZero();
572  return;
573  }
574 
575  /* Let us use the following lemma:
576  *
577  * Lemma: If the matrix A has the LU decomposition PAQ = LU,
578  * then Ker A = Q(Ker U).
579  *
580  * Proof: trivial: just keep in mind that P, Q, L are invertible.
581  */
582 
583  /* Thus, all we need to do is to compute Ker U, and then apply Q.
584  *
585  * U is upper triangular, with eigenvalues sorted so that any zeros appear at the end.
586  * Thus, the diagonal of U ends with exactly
587  * dimKer zero's. Let us use that to construct dimKer linearly
588  * independent vectors in Ker U.
589  */
590 
591  Matrix<Index, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank());
592  RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold();
593  Index p = 0;
594  for(Index i = 0; i < dec().nonzeroPivots(); ++i)
595  if(abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold)
596  pivots.coeffRef(p++) = i;
597  eigen_internal_assert(p == rank());
598 
599  // we construct a temporaty trapezoid matrix m, by taking the U matrix and
600  // permuting the rows and cols to bring the nonnegligible pivots to the top of
601  // the main diagonal. We need that to be able to apply our triangular solvers.
602  // FIXME when we get triangularView-for-rectangular-matrices, this can be simplified
603  Matrix<typename MatrixType::Scalar, Dynamic, Dynamic, MatrixType::Options,
604  MaxSmallDimAtCompileTime, MatrixType::MaxColsAtCompileTime>
605  m(dec().matrixLU().block(0, 0, rank(), cols));
606  for(Index i = 0; i < rank(); ++i)
607  {
608  if(i) m.row(i).head(i).setZero();
609  m.row(i).tail(cols-i) = dec().matrixLU().row(pivots.coeff(i)).tail(cols-i);
610  }
611  m.block(0, 0, rank(), rank());
612  m.block(0, 0, rank(), rank()).template triangularView<StrictlyLower>().setZero();
613  for(Index i = 0; i < rank(); ++i)
614  m.col(i).swap(m.col(pivots.coeff(i)));
615 
616  // ok, we have our trapezoid matrix, we can apply the triangular solver.
617  // notice that the math behind this suggests that we should apply this to the
618  // negative of the RHS, but for performance we just put the negative sign elsewhere, see below.
619  m.topLeftCorner(rank(), rank())
620  .template triangularView<Upper>().solveInPlace(
621  m.topRightCorner(rank(), dimker)
622  );
623 
624  // now we must undo the column permutation that we had applied!
625  for(Index i = rank()-1; i >= 0; --i)
626  m.col(i).swap(m.col(pivots.coeff(i)));
627 
628  // see the negative sign in the next line, that's what we were talking about above.
629  for(Index i = 0; i < rank(); ++i) dst.row(dec().permutationQ().indices().coeff(i)) = -m.row(i).tail(dimker);
630  for(Index i = rank(); i < cols; ++i) dst.row(dec().permutationQ().indices().coeff(i)).setZero();
631  for(Index k = 0; k < dimker; ++k) dst.coeffRef(dec().permutationQ().indices().coeff(rank()+k), k) = Scalar(1);
632  }
633 };
634 
635 /***** Implementation of image() *****************************************************/
636 
637 template<typename _MatrixType>
638 struct image_retval<FullPivLU<_MatrixType> >
639  : image_retval_base<FullPivLU<_MatrixType> >
640 {
641  EIGEN_MAKE_IMAGE_HELPERS(FullPivLU<_MatrixType>)
642 
643  enum { MaxSmallDimAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(
644  MatrixType::MaxColsAtCompileTime,
645  MatrixType::MaxRowsAtCompileTime)
646  };
647 
648  template<typename Dest> void evalTo(Dest& dst) const
649  {
650  if(rank() == 0)
651  {
652  // The Image is just {0}, so it doesn't have a basis properly speaking, but let's
653  // avoid crashing/asserting as that depends on floating point calculations. Let's
654  // just return a single column vector filled with zeros.
655  dst.setZero();
656  return;
657  }
658 
659  Matrix<Index, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank());
660  RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold();
661  Index p = 0;
662  for(Index i = 0; i < dec().nonzeroPivots(); ++i)
663  if(abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold)
664  pivots.coeffRef(p++) = i;
665  eigen_internal_assert(p == rank());
666 
667  for(Index i = 0; i < rank(); ++i)
668  dst.col(i) = originalMatrix().col(dec().permutationQ().indices().coeff(pivots.coeff(i)));
669  }
670 };
671 
672 /***** Implementation of solve() *****************************************************/
673 
674 template<typename _MatrixType, typename Rhs>
675 struct solve_retval<FullPivLU<_MatrixType>, Rhs>
676  : solve_retval_base<FullPivLU<_MatrixType>, Rhs>
677 {
678  EIGEN_MAKE_SOLVE_HELPERS(FullPivLU<_MatrixType>,Rhs)
679 
680  template<typename Dest> void evalTo(Dest& dst) const
681  {
682  /* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1}.
683  * So we proceed as follows:
684  * Step 1: compute c = P * rhs.
685  * Step 2: replace c by the solution x to Lx = c. Exists because L is invertible.
686  * Step 3: replace c by the solution x to Ux = c. May or may not exist.
687  * Step 4: result = Q * c;
688  */
689 
690  const Index rows = dec().rows(), cols = dec().cols(),
691  nonzero_pivots = dec().nonzeroPivots();
692  eigen_assert(rhs().rows() == rows);
693  const Index smalldim = (std::min)(rows, cols);
694 
695  if(nonzero_pivots == 0)
696  {
697  dst.setZero();
698  return;
699  }
700 
701  typename Rhs::PlainObject c(rhs().rows(), rhs().cols());
702 
703  // Step 1
704  c = dec().permutationP() * rhs();
705 
706  // Step 2
707  dec().matrixLU()
708  .topLeftCorner(smalldim,smalldim)
709  .template triangularView<UnitLower>()
710  .solveInPlace(c.topRows(smalldim));
711  if(rows>cols)
712  {
713  c.bottomRows(rows-cols)
714  -= dec().matrixLU().bottomRows(rows-cols)
715  * c.topRows(cols);
716  }
717 
718  // Step 3
719  dec().matrixLU()
720  .topLeftCorner(nonzero_pivots, nonzero_pivots)
721  .template triangularView<Upper>()
722  .solveInPlace(c.topRows(nonzero_pivots));
723 
724  // Step 4
725  for(Index i = 0; i < nonzero_pivots; ++i)
726  dst.row(dec().permutationQ().indices().coeff(i)) = c.row(i);
727  for(Index i = nonzero_pivots; i < dec().matrixLU().cols(); ++i)
728  dst.row(dec().permutationQ().indices().coeff(i)).setZero();
729  }
730 };
731 
732 } // end namespace internal
733 
734 /******* MatrixBase methods *****************************************************************/
735 
742 template<typename Derived>
743 inline const FullPivLU<typename MatrixBase<Derived>::PlainObject>
745 {
746  return FullPivLU<PlainObject>(eval());
747 }
748 
749 } // end namespace Eigen
750 
751 #endif // EIGEN_LU_H