Eigen  3.2.91
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 
15 namespace internal {
16 template<typename _MatrixType> struct traits<FullPivLU<_MatrixType> >
17  : traits<_MatrixType>
18 {
19  enum { Flags = 0 };
20 };
21 
22 } // end namespace internal
23 
55 template<typename _MatrixType> class FullPivLU
56 {
57  public:
58  typedef _MatrixType MatrixType;
59  enum {
60  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
61  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
62  Options = MatrixType::Options,
63  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
64  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
65  };
66  typedef typename MatrixType::Scalar Scalar;
67  typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
68  typedef typename internal::traits<MatrixType>::StorageKind StorageKind;
69  // FIXME should be int
70  typedef typename MatrixType::StorageIndex StorageIndex;
71  typedef typename internal::plain_row_type<MatrixType, StorageIndex>::type IntRowVectorType;
72  typedef typename internal::plain_col_type<MatrixType, StorageIndex>::type IntColVectorType;
73  typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationQType;
74  typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationPType;
75  typedef typename MatrixType::PlainObject PlainObject;
76 
83  FullPivLU();
84 
91  FullPivLU(Index rows, Index cols);
92 
98  explicit FullPivLU(const MatrixType& matrix);
99 
107  FullPivLU& compute(const MatrixType& matrix);
108 
115  inline const MatrixType& matrixLU() const
116  {
117  eigen_assert(m_isInitialized && "LU is not initialized.");
118  return m_lu;
119  }
120 
128  inline Index nonzeroPivots() const
129  {
130  eigen_assert(m_isInitialized && "LU is not initialized.");
131  return m_nonzero_pivots;
132  }
133 
137  RealScalar maxPivot() const { return m_maxpivot; }
138 
143  inline const PermutationPType& permutationP() const
144  {
145  eigen_assert(m_isInitialized && "LU is not initialized.");
146  return m_p;
147  }
148 
153  inline const PermutationQType& permutationQ() const
154  {
155  eigen_assert(m_isInitialized && "LU is not initialized.");
156  return m_q;
157  }
158 
173  inline const internal::kernel_retval<FullPivLU> kernel() const
174  {
175  eigen_assert(m_isInitialized && "LU is not initialized.");
176  return internal::kernel_retval<FullPivLU>(*this);
177  }
178 
198  inline const internal::image_retval<FullPivLU>
199  image(const MatrixType& originalMatrix) const
200  {
201  eigen_assert(m_isInitialized && "LU is not initialized.");
202  return internal::image_retval<FullPivLU>(*this, originalMatrix);
203  }
204 
224  template<typename Rhs>
225  inline const Solve<FullPivLU, Rhs>
226  solve(const MatrixBase<Rhs>& b) const
227  {
228  eigen_assert(m_isInitialized && "LU is not initialized.");
229  return Solve<FullPivLU, Rhs>(*this, b.derived());
230  }
231 
247  typename internal::traits<MatrixType>::Scalar determinant() const;
248 
266  FullPivLU& setThreshold(const RealScalar& threshold)
267  {
268  m_usePrescribedThreshold = true;
269  m_prescribedThreshold = threshold;
270  return *this;
271  }
272 
282  {
283  m_usePrescribedThreshold = false;
284  return *this;
285  }
286 
291  RealScalar threshold() const
292  {
293  eigen_assert(m_isInitialized || m_usePrescribedThreshold);
294  return m_usePrescribedThreshold ? m_prescribedThreshold
295  // this formula comes from experimenting (see "LU precision tuning" thread on the list)
296  // and turns out to be identical to Higham's formula used already in LDLt.
297  : NumTraits<Scalar>::epsilon() * m_lu.diagonalSize();
298  }
299 
306  inline Index rank() const
307  {
308  using std::abs;
309  eigen_assert(m_isInitialized && "LU is not initialized.");
310  RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold();
311  Index result = 0;
312  for(Index i = 0; i < m_nonzero_pivots; ++i)
313  result += (abs(m_lu.coeff(i,i)) > premultiplied_threshold);
314  return result;
315  }
316 
323  inline Index dimensionOfKernel() const
324  {
325  eigen_assert(m_isInitialized && "LU is not initialized.");
326  return cols() - rank();
327  }
328 
336  inline bool isInjective() const
337  {
338  eigen_assert(m_isInitialized && "LU is not initialized.");
339  return rank() == cols();
340  }
341 
349  inline bool isSurjective() const
350  {
351  eigen_assert(m_isInitialized && "LU is not initialized.");
352  return rank() == rows();
353  }
354 
361  inline bool isInvertible() const
362  {
363  eigen_assert(m_isInitialized && "LU is not initialized.");
364  return isInjective() && (m_lu.rows() == m_lu.cols());
365  }
366 
374  inline const Inverse<FullPivLU> inverse() const
375  {
376  eigen_assert(m_isInitialized && "LU is not initialized.");
377  eigen_assert(m_lu.rows() == m_lu.cols() && "You can't take the inverse of a non-square matrix!");
378  return Inverse<FullPivLU>(*this);
379  }
380 
381  MatrixType reconstructedMatrix() const;
382 
383  inline Index rows() const { return m_lu.rows(); }
384  inline Index cols() const { return m_lu.cols(); }
385 
386  #ifndef EIGEN_PARSED_BY_DOXYGEN
387  template<typename RhsType, typename DstType>
388  EIGEN_DEVICE_FUNC
389  void _solve_impl(const RhsType &rhs, DstType &dst) const;
390  #endif
391 
392  protected:
393 
394  static void check_template_parameters()
395  {
396  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
397  }
398 
399  MatrixType m_lu;
400  PermutationPType m_p;
401  PermutationQType m_q;
402  IntColVectorType m_rowsTranspositions;
403  IntRowVectorType m_colsTranspositions;
404  Index m_det_pq, m_nonzero_pivots;
405  RealScalar m_maxpivot, m_prescribedThreshold;
406  bool m_isInitialized, m_usePrescribedThreshold;
407 };
408 
409 template<typename MatrixType>
411  : m_isInitialized(false), m_usePrescribedThreshold(false)
412 {
413 }
414 
415 template<typename MatrixType>
416 FullPivLU<MatrixType>::FullPivLU(Index rows, Index cols)
417  : m_lu(rows, cols),
418  m_p(rows),
419  m_q(cols),
420  m_rowsTranspositions(rows),
421  m_colsTranspositions(cols),
422  m_isInitialized(false),
423  m_usePrescribedThreshold(false)
424 {
425 }
426 
427 template<typename MatrixType>
428 FullPivLU<MatrixType>::FullPivLU(const MatrixType& matrix)
429  : m_lu(matrix.rows(), matrix.cols()),
430  m_p(matrix.rows()),
431  m_q(matrix.cols()),
432  m_rowsTranspositions(matrix.rows()),
433  m_colsTranspositions(matrix.cols()),
434  m_isInitialized(false),
435  m_usePrescribedThreshold(false)
436 {
437  compute(matrix);
438 }
439 
440 template<typename MatrixType>
442 {
443  check_template_parameters();
444 
445  // the permutations are stored as int indices, so just to be sure:
446  eigen_assert(matrix.rows()<=NumTraits<int>::highest() && matrix.cols()<=NumTraits<int>::highest());
447 
448  m_isInitialized = true;
449  m_lu = matrix;
450 
451  const Index size = matrix.diagonalSize();
452  const Index rows = matrix.rows();
453  const Index cols = matrix.cols();
454 
455  // will store the transpositions, before we accumulate them at the end.
456  // can't accumulate on-the-fly because that will be done in reverse order for the rows.
457  m_rowsTranspositions.resize(matrix.rows());
458  m_colsTranspositions.resize(matrix.cols());
459  Index number_of_transpositions = 0; // number of NONTRIVIAL transpositions, i.e. m_rowsTranspositions[i]!=i
460 
461  m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
462  m_maxpivot = RealScalar(0);
463 
464  for(Index k = 0; k < size; ++k)
465  {
466  // First, we need to find the pivot.
467 
468  // biggest coefficient in the remaining bottom-right corner (starting at row k, col k)
469  Index row_of_biggest_in_corner, col_of_biggest_in_corner;
470  typedef internal::scalar_score_coeff_op<Scalar> Scoring;
471  typedef typename Scoring::result_type Score;
472  Score biggest_in_corner;
473  biggest_in_corner = m_lu.bottomRightCorner(rows-k, cols-k)
474  .unaryExpr(Scoring())
475  .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
476  row_of_biggest_in_corner += k; // correct the values! since they were computed in the corner,
477  col_of_biggest_in_corner += k; // need to add k to them.
478 
479  if(biggest_in_corner==Score(0))
480  {
481  // before exiting, make sure to initialize the still uninitialized transpositions
482  // in a sane state without destroying what we already have.
483  m_nonzero_pivots = k;
484  for(Index i = k; i < size; ++i)
485  {
486  m_rowsTranspositions.coeffRef(i) = i;
487  m_colsTranspositions.coeffRef(i) = i;
488  }
489  break;
490  }
491 
492  RealScalar abs_pivot = internal::abs_knowing_score<Scalar>()(m_lu(row_of_biggest_in_corner, col_of_biggest_in_corner), biggest_in_corner);
493  if(abs_pivot > m_maxpivot) m_maxpivot = abs_pivot;
494 
495  // Now that we've found the pivot, we need to apply the row/col swaps to
496  // bring it to the location (k,k).
497 
498  m_rowsTranspositions.coeffRef(k) = row_of_biggest_in_corner;
499  m_colsTranspositions.coeffRef(k) = col_of_biggest_in_corner;
500  if(k != row_of_biggest_in_corner) {
501  m_lu.row(k).swap(m_lu.row(row_of_biggest_in_corner));
502  ++number_of_transpositions;
503  }
504  if(k != col_of_biggest_in_corner) {
505  m_lu.col(k).swap(m_lu.col(col_of_biggest_in_corner));
506  ++number_of_transpositions;
507  }
508 
509  // Now that the pivot is at the right location, we update the remaining
510  // bottom-right corner by Gaussian elimination.
511 
512  if(k<rows-1)
513  m_lu.col(k).tail(rows-k-1) /= m_lu.coeff(k,k);
514  if(k<size-1)
515  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);
516  }
517 
518  // the main loop is over, we still have to accumulate the transpositions to find the
519  // permutations P and Q
520 
521  m_p.setIdentity(rows);
522  for(Index k = size-1; k >= 0; --k)
523  m_p.applyTranspositionOnTheRight(k, m_rowsTranspositions.coeff(k));
524 
525  m_q.setIdentity(cols);
526  for(Index k = 0; k < size; ++k)
527  m_q.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k));
528 
529  m_det_pq = (number_of_transpositions%2) ? -1 : 1;
530  return *this;
531 }
532 
533 template<typename MatrixType>
534 typename internal::traits<MatrixType>::Scalar FullPivLU<MatrixType>::determinant() const
535 {
536  eigen_assert(m_isInitialized && "LU is not initialized.");
537  eigen_assert(m_lu.rows() == m_lu.cols() && "You can't take the determinant of a non-square matrix!");
538  return Scalar(m_det_pq) * Scalar(m_lu.diagonal().prod());
539 }
540 
544 template<typename MatrixType>
546 {
547  eigen_assert(m_isInitialized && "LU is not initialized.");
548  const Index smalldim = (std::min)(m_lu.rows(), m_lu.cols());
549  // LU
550  MatrixType res(m_lu.rows(),m_lu.cols());
551  // FIXME the .toDenseMatrix() should not be needed...
552  res = m_lu.leftCols(smalldim)
553  .template triangularView<UnitLower>().toDenseMatrix()
554  * m_lu.topRows(smalldim)
555  .template triangularView<Upper>().toDenseMatrix();
556 
557  // P^{-1}(LU)
558  res = m_p.inverse() * res;
559 
560  // (P^{-1}LU)Q^{-1}
561  res = res * m_q.inverse();
562 
563  return res;
564 }
565 
566 /********* Implementation of kernel() **************************************************/
567 
568 namespace internal {
569 template<typename _MatrixType>
570 struct kernel_retval<FullPivLU<_MatrixType> >
571  : kernel_retval_base<FullPivLU<_MatrixType> >
572 {
573  EIGEN_MAKE_KERNEL_HELPERS(FullPivLU<_MatrixType>)
574 
575  enum { MaxSmallDimAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(
576  MatrixType::MaxColsAtCompileTime,
577  MatrixType::MaxRowsAtCompileTime)
578  };
579 
580  template<typename Dest> void evalTo(Dest& dst) const
581  {
582  using std::abs;
583  const Index cols = dec().matrixLU().cols(), dimker = cols - rank();
584  if(dimker == 0)
585  {
586  // The Kernel is just {0}, so it doesn't have a basis properly speaking, but let's
587  // avoid crashing/asserting as that depends on floating point calculations. Let's
588  // just return a single column vector filled with zeros.
589  dst.setZero();
590  return;
591  }
592 
593  /* Let us use the following lemma:
594  *
595  * Lemma: If the matrix A has the LU decomposition PAQ = LU,
596  * then Ker A = Q(Ker U).
597  *
598  * Proof: trivial: just keep in mind that P, Q, L are invertible.
599  */
600 
601  /* Thus, all we need to do is to compute Ker U, and then apply Q.
602  *
603  * U is upper triangular, with eigenvalues sorted so that any zeros appear at the end.
604  * Thus, the diagonal of U ends with exactly
605  * dimKer zero's. Let us use that to construct dimKer linearly
606  * independent vectors in Ker U.
607  */
608 
609  Matrix<Index, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank());
610  RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold();
611  Index p = 0;
612  for(Index i = 0; i < dec().nonzeroPivots(); ++i)
613  if(abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold)
614  pivots.coeffRef(p++) = i;
615  eigen_internal_assert(p == rank());
616 
617  // we construct a temporaty trapezoid matrix m, by taking the U matrix and
618  // permuting the rows and cols to bring the nonnegligible pivots to the top of
619  // the main diagonal. We need that to be able to apply our triangular solvers.
620  // FIXME when we get triangularView-for-rectangular-matrices, this can be simplified
621  Matrix<typename MatrixType::Scalar, Dynamic, Dynamic, MatrixType::Options,
622  MaxSmallDimAtCompileTime, MatrixType::MaxColsAtCompileTime>
623  m(dec().matrixLU().block(0, 0, rank(), cols));
624  for(Index i = 0; i < rank(); ++i)
625  {
626  if(i) m.row(i).head(i).setZero();
627  m.row(i).tail(cols-i) = dec().matrixLU().row(pivots.coeff(i)).tail(cols-i);
628  }
629  m.block(0, 0, rank(), rank());
630  m.block(0, 0, rank(), rank()).template triangularView<StrictlyLower>().setZero();
631  for(Index i = 0; i < rank(); ++i)
632  m.col(i).swap(m.col(pivots.coeff(i)));
633 
634  // ok, we have our trapezoid matrix, we can apply the triangular solver.
635  // notice that the math behind this suggests that we should apply this to the
636  // negative of the RHS, but for performance we just put the negative sign elsewhere, see below.
637  m.topLeftCorner(rank(), rank())
638  .template triangularView<Upper>().solveInPlace(
639  m.topRightCorner(rank(), dimker)
640  );
641 
642  // now we must undo the column permutation that we had applied!
643  for(Index i = rank()-1; i >= 0; --i)
644  m.col(i).swap(m.col(pivots.coeff(i)));
645 
646  // see the negative sign in the next line, that's what we were talking about above.
647  for(Index i = 0; i < rank(); ++i) dst.row(dec().permutationQ().indices().coeff(i)) = -m.row(i).tail(dimker);
648  for(Index i = rank(); i < cols; ++i) dst.row(dec().permutationQ().indices().coeff(i)).setZero();
649  for(Index k = 0; k < dimker; ++k) dst.coeffRef(dec().permutationQ().indices().coeff(rank()+k), k) = Scalar(1);
650  }
651 };
652 
653 /***** Implementation of image() *****************************************************/
654 
655 template<typename _MatrixType>
656 struct image_retval<FullPivLU<_MatrixType> >
657  : image_retval_base<FullPivLU<_MatrixType> >
658 {
659  EIGEN_MAKE_IMAGE_HELPERS(FullPivLU<_MatrixType>)
660 
661  enum { MaxSmallDimAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(
662  MatrixType::MaxColsAtCompileTime,
663  MatrixType::MaxRowsAtCompileTime)
664  };
665 
666  template<typename Dest> void evalTo(Dest& dst) const
667  {
668  using std::abs;
669  if(rank() == 0)
670  {
671  // The Image is just {0}, so it doesn't have a basis properly speaking, but let's
672  // avoid crashing/asserting as that depends on floating point calculations. Let's
673  // just return a single column vector filled with zeros.
674  dst.setZero();
675  return;
676  }
677 
678  Matrix<Index, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank());
679  RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold();
680  Index p = 0;
681  for(Index i = 0; i < dec().nonzeroPivots(); ++i)
682  if(abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold)
683  pivots.coeffRef(p++) = i;
684  eigen_internal_assert(p == rank());
685 
686  for(Index i = 0; i < rank(); ++i)
687  dst.col(i) = originalMatrix().col(dec().permutationQ().indices().coeff(pivots.coeff(i)));
688  }
689 };
690 
691 /***** Implementation of solve() *****************************************************/
692 
693 } // end namespace internal
694 
695 #ifndef EIGEN_PARSED_BY_DOXYGEN
696 template<typename _MatrixType>
697 template<typename RhsType, typename DstType>
698 void FullPivLU<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const
699 {
700  /* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1}.
701  * So we proceed as follows:
702  * Step 1: compute c = P * rhs.
703  * Step 2: replace c by the solution x to Lx = c. Exists because L is invertible.
704  * Step 3: replace c by the solution x to Ux = c. May or may not exist.
705  * Step 4: result = Q * c;
706  */
707 
708  const Index rows = this->rows(),
709  cols = this->cols(),
710  nonzero_pivots = this->nonzeroPivots();
711  eigen_assert(rhs.rows() == rows);
712  const Index smalldim = (std::min)(rows, cols);
713 
714  if(nonzero_pivots == 0)
715  {
716  dst.setZero();
717  return;
718  }
719 
720  typename RhsType::PlainObject c(rhs.rows(), rhs.cols());
721 
722  // Step 1
723  c = permutationP() * rhs;
724 
725  // Step 2
726  m_lu.topLeftCorner(smalldim,smalldim)
727  .template triangularView<UnitLower>()
728  .solveInPlace(c.topRows(smalldim));
729  if(rows>cols)
730  c.bottomRows(rows-cols) -= m_lu.bottomRows(rows-cols) * c.topRows(cols);
731 
732  // Step 3
733  m_lu.topLeftCorner(nonzero_pivots, nonzero_pivots)
734  .template triangularView<Upper>()
735  .solveInPlace(c.topRows(nonzero_pivots));
736 
737  // Step 4
738  for(Index i = 0; i < nonzero_pivots; ++i)
739  dst.row(permutationQ().indices().coeff(i)) = c.row(i);
740  for(Index i = nonzero_pivots; i < m_lu.cols(); ++i)
741  dst.row(permutationQ().indices().coeff(i)).setZero();
742 }
743 #endif
744 
745 namespace internal {
746 
747 
748 /***** Implementation of inverse() *****************************************************/
749 template<typename DstXprType, typename MatrixType, typename Scalar>
750 struct Assignment<DstXprType, Inverse<FullPivLU<MatrixType> >, internal::assign_op<Scalar>, Dense2Dense, Scalar>
751 {
752  typedef FullPivLU<MatrixType> LuType;
753  typedef Inverse<LuType> SrcXprType;
754  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar> &)
755  {
756  dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
757  }
758 };
759 } // end namespace internal
760 
761 /******* MatrixBase methods *****************************************************************/
762 
769 #ifndef __CUDACC__
770 template<typename Derived>
771 inline const FullPivLU<typename MatrixBase<Derived>::PlainObject>
773 {
774  return FullPivLU<PlainObject>(eval());
775 }
776 #endif
777 
778 } // end namespace Eigen
779 
780 #endif // EIGEN_LU_H
bool isInjective() const
Definition: FullPivLU.h:336
const internal::image_retval< FullPivLU > image(const MatrixType &originalMatrix) const
Definition: FullPivLU.h:199
FullPivLU & setThreshold(const RealScalar &threshold)
Definition: FullPivLU.h:266
internal::traits< MatrixType >::Scalar determinant() const
Definition: FullPivLU.h:534
Definition: LDLT.h:16
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:107
RealScalar threshold() const
Definition: FullPivLU.h:291
Index nonzeroPivots() const
Definition: FullPivLU.h:128
Expression of the inverse of another expression.
Definition: Inverse.h:45
bool isSurjective() const
Definition: FullPivLU.h:349
FullPivLU & setThreshold(Default_t)
Definition: FullPivLU.h:281
const PermutationQType & permutationQ() const
Definition: FullPivLU.h:153
FullPivLU & compute(const MatrixType &matrix)
Definition: FullPivLU.h:441
MatrixType reconstructedMatrix() const
Definition: FullPivLU.h:545
const PermutationPType & permutationP() const
Definition: FullPivLU.h:143
Definition: Eigen_Colamd.h:54
LU decomposition of a matrix with complete pivoting, and related features.
Definition: ForwardDeclarations.h:246
bool isInvertible() const
Definition: FullPivLU.h:361
const FullPivLU< PlainObject > fullPivLu() const
Definition: FullPivLU.h:772
FullPivLU()
Default Constructor.
Definition: FullPivLU.h:410
Index rank() const
Definition: FullPivLU.h:306
const MatrixType & matrixLU() const
Definition: FullPivLU.h:115
Pseudo expression representing a solving operation.
Definition: Solve.h:63
const Inverse< FullPivLU > inverse() const
Definition: FullPivLU.h:374
Index dimensionOfKernel() const
Definition: FullPivLU.h:323
RealScalar maxPivot() const
Definition: FullPivLU.h:137
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
const Solve< FullPivLU, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: FullPivLU.h:226
const internal::kernel_retval< FullPivLU > kernel() const
Definition: FullPivLU.h:173