GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-07-11 07:31:48 +08:00
Code Issues Packages Projects Releases Wiki Activity

bug #941: fix accuracy issue in ColPivHouseholderQR, do not stop decomposition on a small pivot

Browse Source 
(grafted from f1092d2f736bbe541b3d8b5cca893f668f9c6b5f
)
This commit is contained in:
Gael Guennebaud 2015-01-30 19:04:04 +01:00
parent 8296c4aaed
commit f89ba2a58b
1 changed files with 8 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
Eigen/src/QR/ColPivHouseholderQR.h
View File

@ -463,20 +463,10 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const
// we store that back into our table: it can't hurt to correct our table. // we store that back into our table: it can't hurt to correct our table.
m_colSqNorms.coeffRef(biggest_col_index) = biggest_col_sq_norm; m_colSqNorms.coeffRef(biggest_col_index) = biggest_col_sq_norm;
// if the current biggest column is smaller than epsilon times the initial biggest column, // Track the number of meaningful pivots but do not stop the decomposition to make
// terminate to avoid generating nan/inf values. // sure that the initial matrix is properly reproduced. See bug 941.
// Note that here, if we test instead for "biggest == 0", we get a failure every 1000 (or so) if(m_nonzero_pivots==size && biggest_col_sq_norm < threshold_helper * RealScalar(rows-k))
// repetitions of the unit test, with the result of solve() filled with large values of the order
// of 1/(size*epsilon).
if(biggest_col_sq_norm < threshold_helper * RealScalar(rows-k))
{
m_nonzero_pivots = k; m_nonzero_pivots = k;
m_hCoeffs.tail(size-k).setZero();
m_qr.bottomRightCorner(rows-k,cols-k)
.template triangularView<StrictlyLower>()
.setZero();
break;
}
// apply the transposition to the columns // apply the transposition to the columns
m_colsTranspositions.coeffRef(k) = biggest_col_index; m_colsTranspositions.coeffRef(k) = biggest_col_index;
@ -505,7 +495,7 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const
} }
m_colsPermutation.setIdentity(PermIndexType(cols)); m_colsPermutation.setIdentity(PermIndexType(cols));
for(PermIndexType k = 0; k < m_nonzero_pivots; ++k) for(PermIndexType k = 0; k < size/*m_nonzero_pivots*/; ++k)
m_colsPermutation.applyTranspositionOnTheRight(k, PermIndexType(m_colsTranspositions.coeff(k))); m_colsPermutation.applyTranspositionOnTheRight(k, PermIndexType(m_colsTranspositions.coeff(k)));
m_det_pq = (number_of_transpositions%2) ? -1 : 1; m_det_pq = (number_of_transpositions%2) ? -1 : 1;
@ -555,13 +545,15 @@ struct solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs>
} // end namespace internal } // end namespace internal
/** \returns the matrix Q as a sequence of householder transformations */ /** \returns the matrix Q as a sequence of householder transformations.
* You can extract the meaningful part only by using:
* \code qr.householderQ().setLength(qr.nonzeroPivots()) */
template<typename MatrixType> template<typename MatrixType>
typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType> typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType>
::householderQ() const ::householderQ() const
{ {
eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate()).setLength(m_nonzero_pivots); return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
} }
/** \return the column-pivoting Householder QR decomposition of \c *this. /** \return the column-pivoting Householder QR decomposition of \c *this.