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

bug #1214: consider denormals as zero in D&C SVD. This also workaround infinite binary search when compiling with ICC's unsafe optimizations.

Browse Source
This commit is contained in:
Gael Guennebaud 2016-05-03 23:15:29 +02:00
parent 2c5568a757
commit e2ca478485
1 changed files with 24 additions and 13 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

37
Eigen/src/SVD/BDCSVD.h
View File

@ -11,7 +11,7 @@
// Copyright (C) 2013 Jean Ceccato <jean.ceccato@ensimag.fr> // Copyright (C) 2013 Jean Ceccato <jean.ceccato@ensimag.fr>
// Copyright (C) 2013 Pierre Zoppitelli <pierre.zoppitelli@ensimag.fr> // Copyright (C) 2013 Pierre Zoppitelli <pierre.zoppitelli@ensimag.fr>
// Copyright (C) 2013 Jitse Niesen <jitse@maths.leeds.ac.uk> // Copyright (C) 2013 Jitse Niesen <jitse@maths.leeds.ac.uk>
// Copyright (C) 2014 Gael Guennebaud <gael.guennebaud@inria.fr> // Copyright (C) 2014-2016 Gael Guennebaud <gael.guennebaud@inria.fr>
// //
// Source Code Form is subject to the terms of the Mozilla // Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed // Public License v. 2.0. If a copy of the MPL was not distributed
@ -21,6 +21,7 @@
#define EIGEN_BDCSVD_H #define EIGEN_BDCSVD_H
// #define EIGEN_BDCSVD_DEBUG_VERBOSE // #define EIGEN_BDCSVD_DEBUG_VERBOSE
// #define EIGEN_BDCSVD_SANITY_CHECKS // #define EIGEN_BDCSVD_SANITY_CHECKS
namespace Eigen { namespace Eigen {
#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
@ -228,6 +229,8 @@ BDCSVD<MatrixType>& BDCSVD<MatrixType>::compute(const MatrixType& matrix, unsign
#endif #endif
allocate(matrix.rows(), matrix.cols(), computationOptions); allocate(matrix.rows(), matrix.cols(), computationOptions);
using std::abs; using std::abs;
const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)();
//**** step -1 - If the problem is too small, directly falls back to JacobiSVD and return //**** step -1 - If the problem is too small, directly falls back to JacobiSVD and return
if(matrix.cols() < m_algoswap) if(matrix.cols() < m_algoswap)
@ -266,7 +269,7 @@ BDCSVD<MatrixType>& BDCSVD<MatrixType>::compute(const MatrixType& matrix, unsign
{ {
RealScalar a = abs(m_computed.coeff(i, i)); RealScalar a = abs(m_computed.coeff(i, i));
m_singularValues.coeffRef(i) = a * scale; m_singularValues.coeffRef(i) = a * scale;
if (a == 0) if (a<considerZero)
{ {
m_nonzeroSingularValues = i; m_nonzeroSingularValues = i;
m_singularValues.tail(m_diagSize - i - 1).setZero(); m_singularValues.tail(m_diagSize - i - 1).setZero();
@ -380,6 +383,7 @@ void BDCSVD<MatrixType>::divide (Index firstCol, Index lastCol, Index firstRowW,
using std::abs; using std::abs;
const Index n = lastCol - firstCol + 1; const Index n = lastCol - firstCol + 1;
const Index k = n/2; const Index k = n/2;
const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)();
RealScalar alphaK; RealScalar alphaK;
RealScalar betaK; RealScalar betaK;
RealScalar r0; RealScalar r0;
@ -434,7 +438,7 @@ void BDCSVD<MatrixType>::divide (Index firstCol, Index lastCol, Index firstRowW,
f = m_naiveU.row(0).segment(firstCol + k + 1, n - k - 1); f = m_naiveU.row(0).segment(firstCol + k + 1, n - k - 1);
} }
if (m_compV) m_naiveV(firstRowW+k, firstColW) = 1; if (m_compV) m_naiveV(firstRowW+k, firstColW) = 1;
if (r0 == 0) if (r0<considerZero)
{ {
c0 = 1; c0 = 1;
s0 = 0; s0 = 0;
@ -553,6 +557,8 @@ void BDCSVD<MatrixType>::divide (Index firstCol, Index lastCol, Index firstRowW,
template <typename MatrixType> template <typename MatrixType>
void BDCSVD<MatrixType>::computeSVDofM(Index firstCol, Index n, MatrixXr& U, VectorType& singVals, MatrixXr& V) void BDCSVD<MatrixType>::computeSVDofM(Index firstCol, Index n, MatrixXr& U, VectorType& singVals, MatrixXr& V)
{ {
const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)();
using std::abs;
ArrayRef col0 = m_computed.col(firstCol).segment(firstCol, n); ArrayRef col0 = m_computed.col(firstCol).segment(firstCol, n);
m_workspace.head(n) = m_computed.block(firstCol, firstCol, n, n).diagonal(); m_workspace.head(n) = m_computed.block(firstCol, firstCol, n, n).diagonal();
ArrayRef diag = m_workspace.head(n); ArrayRef diag = m_workspace.head(n);
@ -575,7 +581,7 @@ void BDCSVD<MatrixType>::computeSVDofM(Index firstCol, Index n, MatrixXr& U, Vec
while(actual_n>1 && diag(actual_n-1)==0) --actual_n; while(actual_n>1 && diag(actual_n-1)==0) --actual_n;
Index m = 0; // size of the deflated problem Index m = 0; // size of the deflated problem
for(Index k=0;k<actual_n;++k) for(Index k=0;k<actual_n;++k)
if(col0(k)!=0) if(abs(col0(k))>considerZero)
m_workspaceI(m++) = k; m_workspaceI(m++) = k;
Map<ArrayXi> perm(m_workspaceI.data(),m); Map<ArrayXi> perm(m_workspaceI.data(),m);
@ -600,7 +606,7 @@ void BDCSVD<MatrixType>::computeSVDofM(Index firstCol, Index n, MatrixXr& U, Vec
{ {
Index actual_n = n; Index actual_n = n;
while(actual_n>1 && col0(actual_n-1)==0) --actual_n; while(actual_n>1 && abs(col0(actual_n-1))<considerZero) --actual_n;
std::cout << "\n\n mus: " << mus.head(actual_n).transpose() << "\n\n"; std::cout << "\n\n mus: " << mus.head(actual_n).transpose() << "\n\n";
std::cout << " check1 (expect0) : " << ((singVals.array()-(shifts+mus)) / singVals.array()).head(actual_n).transpose() << "\n\n"; std::cout << " check1 (expect0) : " << ((singVals.array()-(shifts+mus)) / singVals.array()).head(actual_n).transpose() << "\n\n";
std::cout << " check2 (>0) : " << ((singVals.array()-diag) / singVals.array()).head(actual_n).transpose() << "\n\n"; std::cout << " check2 (>0) : " << ((singVals.array()-diag) / singVals.array()).head(actual_n).transpose() << "\n\n";
@ -680,6 +686,7 @@ typename BDCSVD<MatrixType>::RealScalar BDCSVD<MatrixType>::secularEq(RealScalar
res += numext::abs2(col0(j)) / ((diagShifted(j) - mu) * (diag(j) + shift + mu)); res += numext::abs2(col0(j)) / ((diagShifted(j) - mu) * (diag(j) + shift + mu));
} }
return res; return res;
} }
template <typename MatrixType> template <typename MatrixType>
@ -798,7 +805,7 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
RealScalar leftShifted, rightShifted; RealScalar leftShifted, rightShifted;
if (shift == left) if (shift == left)
{ {
leftShifted = RealScalar(1)/NumTraits<RealScalar>::highest(); leftShifted = (std::numeric_limits<RealScalar>::min)();
// I don't understand why the case k==0 would be special there: // I don't understand why the case k==0 would be special there:
// if (k == 0) rightShifted = right - left; else // if (k == 0) rightShifted = right - left; else
rightShifted = (k==actual_n-1) ? right : ((right - left) * 0.6); // theoretically we can take 0.5, but let's be safe rightShifted = (k==actual_n-1) ? right : ((right - left) * 0.6); // theoretically we can take 0.5, but let's be safe
@ -806,7 +813,7 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
else else
{ {
leftShifted = -(right - left) * 0.6; leftShifted = -(right - left) * 0.6;
rightShifted = -RealScalar(1)/NumTraits<RealScalar>::highest(); rightShifted = -(std::numeric_limits<RealScalar>::min)();
} }
RealScalar fLeft = secularEq(leftShifted, col0, diag, perm, diagShifted, shift); RealScalar fLeft = secularEq(leftShifted, col0, diag, perm, diagShifted, shift);
@ -817,7 +824,10 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
if(!(fLeft * fRight<0)) if(!(fLeft * fRight<0))
{
std::cout << "fLeft: " << leftShifted << " - " << diagShifted.head(10).transpose() << "\n ; " << bool(left==shift) << " " << (left-shift) << "\n";
std::cout << k << " : " << fLeft << " * " << fRight << " == " << fLeft * fRight << " ; " << left << " - " << right << " -> " << leftShifted << " " << rightShifted << " shift=" << shift << "\n"; std::cout << k << " : " << fLeft << " * " << fRight << " == " << fLeft * fRight << " ; " << left << " - " << right << " -> " << leftShifted << " " << rightShifted << " shift=" << shift << "\n";
}
#endif #endif
eigen_internal_assert(fLeft * fRight < 0); eigen_internal_assert(fLeft * fRight < 0);
@ -1028,8 +1038,9 @@ void BDCSVD<MatrixType>::deflation(Index firstCol, Index lastCol, Index k, Index
Diagonal<MatrixXr> fulldiag(m_computed); Diagonal<MatrixXr> fulldiag(m_computed);
VectorBlock<Diagonal<MatrixXr>,Dynamic> diag(fulldiag, firstCol+shift, length); VectorBlock<Diagonal<MatrixXr>,Dynamic> diag(fulldiag, firstCol+shift, length);
const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)();
RealScalar maxDiag = diag.tail((std::max)(Index(1),length-1)).cwiseAbs().maxCoeff(); RealScalar maxDiag = diag.tail((std::max)(Index(1),length-1)).cwiseAbs().maxCoeff();
RealScalar epsilon_strict = NumTraits<RealScalar>::epsilon() * maxDiag; RealScalar epsilon_strict = numext::maxi(considerZero,NumTraits<RealScalar>::epsilon() * maxDiag);
RealScalar epsilon_coarse = 8 * NumTraits<RealScalar>::epsilon() * numext::maxi<RealScalar>(col0.cwiseAbs().maxCoeff(), maxDiag); RealScalar epsilon_coarse = 8 * NumTraits<RealScalar>::epsilon() * numext::maxi<RealScalar>(col0.cwiseAbs().maxCoeff(), maxDiag);
#ifdef EIGEN_BDCSVD_SANITY_CHECKS #ifdef EIGEN_BDCSVD_SANITY_CHECKS
@ -1082,7 +1093,7 @@ void BDCSVD<MatrixType>::deflation(Index firstCol, Index lastCol, Index k, Index
{ {
// Check for total deflation // Check for total deflation
// If we have a total deflation, then we have to consider col0(0)==diag(0) as a singular value during sorting // If we have a total deflation, then we have to consider col0(0)==diag(0) as a singular value during sorting
bool total_deflation = (col0.tail(length-1).array()==RealScalar(0)).all(); bool total_deflation = (col0.tail(length-1).array()<considerZero).all();
// Sort the diagonal entries, since diag(1:k-1) and diag(k:length) are already sorted, let's do a sorted merge. // Sort the diagonal entries, since diag(1:k-1) and diag(k:length) are already sorted, let's do a sorted merge.
// First, compute the respective permutation. // First, compute the respective permutation.
@ -1093,7 +1104,7 @@ void BDCSVD<MatrixType>::deflation(Index firstCol, Index lastCol, Index k, Index
// Move deflated diagonal entries at the end. // Move deflated diagonal entries at the end.
for(Index i=1; i<length; ++i) for(Index i=1; i<length; ++i)
if(diag(i)==0) if(abs(diag(i))<considerZero)
permutation[p++] = i; permutation[p++] = i;
Index i=1, j=k+1; Index i=1, j=k+1;
@ -1112,7 +1123,7 @@ void BDCSVD<MatrixType>::deflation(Index firstCol, Index lastCol, Index k, Index
for(Index i=1; i<length; ++i) for(Index i=1; i<length; ++i)
{ {
Index pi = permutation[i]; Index pi = permutation[i];
if(diag(pi)==0 || diag(0)<diag(pi)) if(abs(diag(pi))<considerZero || diag(0)<diag(pi))
permutation[i-1] = permutation[i]; permutation[i-1] = permutation[i];
else else
{ {
@ -1163,7 +1174,7 @@ void BDCSVD<MatrixType>::deflation(Index firstCol, Index lastCol, Index k, Index
//condition 4.4 //condition 4.4
{ {
Index i = length-1; Index i = length-1;
while(i>0 && (diag(i)==0 || col0(i)==0)) --i; while(i>0 && (abs(diag(i))<considerZero || abs(col0(i))<considerZero)) --i;
for(; i>1;--i) for(; i>1;--i)
if( (diag(i) - diag(i-1)) < NumTraits<RealScalar>::epsilon()*maxDiag ) if( (diag(i) - diag(i-1)) < NumTraits<RealScalar>::epsilon()*maxDiag )
{ {
@ -1177,7 +1188,7 @@ void BDCSVD<MatrixType>::deflation(Index firstCol, Index lastCol, Index k, Index
#ifdef EIGEN_BDCSVD_SANITY_CHECKS #ifdef EIGEN_BDCSVD_SANITY_CHECKS
for(Index j=2;j<length;++j) for(Index j=2;j<length;++j)
assert(diag(j-1)<=diag(j) || diag(j)==0); assert(diag(j-1)<=diag(j) || abs(diag(j))<considerZero);
#endif #endif
#ifdef EIGEN_BDCSVD_SANITY_CHECKS #ifdef EIGEN_BDCSVD_SANITY_CHECKS