mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-03-01 18:26:24 +08:00
Code Issues Packages Projects Releases Wiki Activity

clarifications in LU::solve() and in LU documentation

Browse Source
This commit is contained in:
Benoit Jacob 2009-08-24 00:02:49 -04:00
parent 0926549659
commit b37ab9b324
1 changed files with 14 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
Eigen/src/LU/LU.h
View File

@ -42,11 +42,10 @@
* This decomposition provides the generic approach to solving systems of linear equations, computing * This decomposition provides the generic approach to solving systems of linear equations, computing
* the rank, invertibility, inverse, kernel, and determinant. * the rank, invertibility, inverse, kernel, and determinant.
* *
* This LU decomposition is very stable and well tested with large matrices. Even exact rank computation * This LU decomposition is very stable and well tested with large matrices. However there are use cases where the SVD
* works at sizes larger than 1000x1000. However there are use cases where the SVD decomposition is inherently * decomposition is inherently more stable and/or flexible. For example, when computing the kernel of a matrix,
* more stable when dealing with numerically damaged input. For example, computing the kernel is more stable with * working with the SVD allows to select the smallest singular values of the matrix, something that
* SVD because the SVD can determine which singular values are negligible while LU has to work at the level of matrix * the LU decomposition doesn't see.
* coefficients that are less meaningful in this respect.
* *
* The data of the LU decomposition can be directly accessed through the methods matrixLU(), * The data of the LU decomposition can be directly accessed through the methods matrixLU(),
* permutationP(), permutationQ(). * permutationP(), permutationQ().
@ -326,6 +325,7 @@ template<typename MatrixType> class LU
inline void computeInverse(MatrixType *result) const inline void computeInverse(MatrixType *result) const
{ {
ei_assert(m_originalMatrix != 0 && "LU is not initialized."); ei_assert(m_originalMatrix != 0 && "LU is not initialized.");
ei_assert(m_lu.rows() == m_lu.cols() && "You can't take the inverse of a non-square matrix!");
solve(MatrixType::Identity(m_lu.rows(), m_lu.cols()), result); solve(MatrixType::Identity(m_lu.rows(), m_lu.cols()), result);
} }
@ -456,6 +456,7 @@ template<typename MatrixType>
typename ei_traits<MatrixType>::Scalar LU<MatrixType>::determinant() const typename ei_traits<MatrixType>::Scalar LU<MatrixType>::determinant() const
{ {
ei_assert(m_originalMatrix != 0 && "LU is not initialized."); ei_assert(m_originalMatrix != 0 && "LU is not initialized.");
ei_assert(m_lu.rows() == m_lu.cols() && "You can't take the determinant of a non-square matrix!");
return Scalar(m_det_pq) * m_lu.diagonal().prod(); return Scalar(m_det_pq) * m_lu.diagonal().prod();
} }
@ -533,7 +534,8 @@ bool LU<MatrixType>::solve(
) const ) const
{ {
ei_assert(m_originalMatrix != 0 && "LU is not initialized."); ei_assert(m_originalMatrix != 0 && "LU is not initialized.");
if(m_rank==0) return false;
/* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1}. /* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1}.
* So we proceed as follows: * So we proceed as follows:
* Step 1: compute c = Pb. * Step 1: compute c = Pb.
@ -552,7 +554,8 @@ bool LU<MatrixType>::solve(
for(int i = 0; i < rows; ++i) c.row(m_p.coeff(i)) = b.row(i); for(int i = 0; i < rows; ++i) c.row(m_p.coeff(i)) = b.row(i);
// Step 2 // Step 2
m_lu.corner(Eigen::TopLeft,smalldim,smalldim).template triangularView<UnitLowerTriangular>() m_lu.corner(Eigen::TopLeft,smalldim,smalldim)
.template triangularView<UnitLowerTriangular>()
.solveInPlace(c.corner(Eigen::TopLeft, smalldim, c.cols())); .solveInPlace(c.corner(Eigen::TopLeft, smalldim, c.cols()));
if(rows>cols) if(rows>cols)
{ {
@ -564,11 +567,10 @@ bool LU<MatrixType>::solve(
if(!isSurjective()) if(!isSurjective())
{ {
// is c is in the image of U ? // is c is in the image of U ?
RealScalar biggest_in_c = m_rank>0 ? c.corner(TopLeft, m_rank, c.cols()).cwise().abs().maxCoeff() : 0; RealScalar biggest_in_upper_part_of_c = c.corner(TopLeft, m_rank, c.cols()).cwise().abs().maxCoeff();
for(int col = 0; col < c.cols(); ++col) RealScalar biggest_in_lower_part_of_c = c.corner(BottomLeft, rows-m_rank, c.cols()).cwise().abs().maxCoeff();
for(int row = m_rank; row < c.rows(); ++row) if(!ei_isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision))
if(!ei_isMuchSmallerThan(c.coeff(row,col), biggest_in_c, m_precision)) return false;
return false;
} }
m_lu.corner(TopLeft, m_rank, m_rank) m_lu.corner(TopLeft, m_rank, m_rank)
.template triangularView<UpperTriangular>() .template triangularView<UpperTriangular>()