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


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Mon Jun 10 2019 12:34:36