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

bug #1695: fix a numerical robustness issue. Computing the secular equation at the middle range without a shift might give a wrong sign.

Browse Source
This commit is contained in:
Gael Guennebaud 2019-03-27 20:16:58 +01:00
parent 8de66719f9
commit 45e65fbb77
1 changed files with 54 additions and 24 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -785,6 +785,21 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
Map<ArrayXr> diagShifted(m_workspace.data()+4*n, n); Map<ArrayXr> diagShifted(m_workspace.data()+4*n, n);
diagShifted = diag - shift; diagShifted = diag - shift;
if(k!=actual_n-1)
{
// check that after the shift, f(mid) is still negative:
RealScalar midShifted = (right - left) / RealScalar(2);
if(shift==right)
midShifted = -midShifted;
RealScalar fMidShifted = secularEq(midShifted, col0, diag, perm, diagShifted, shift);
if(fMidShifted>0)
{
// fMid was erroneous, fix it:
shift = fMidShifted > Literal(0) ? left : right;
diagShifted = diag - shift;
}
}
// initial guess // initial guess
RealScalar muPrev, muCur; RealScalar muPrev, muCur;
if (shift == left) if (shift == left)
@ -822,7 +837,7 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
RealScalar fZero = secularEq(muZero, col0, diag, perm, diagShifted, shift); RealScalar fZero = secularEq(muZero, col0, diag, perm, diagShifted, shift);
#ifdef EIGEN_BDCSVD_SANITY_CHECKS #ifdef EIGEN_BDCSVD_SANITY_CHECKS
assert((std::isfinite)(fZero)); assert((numext::isfinite)(fZero));
#endif #endif
muPrev = muCur; muPrev = muCur;
@ -864,19 +879,20 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
} }
RealScalar fLeft = secularEq(leftShifted, col0, diag, perm, diagShifted, shift); RealScalar fLeft = secularEq(leftShifted, col0, diag, perm, diagShifted, shift);
eigen_internal_assert(fLeft<Literal(0));
#if defined EIGEN_INTERNAL_DEBUGGING || defined EIGEN_BDCSVD_SANITY_CHECKS #if defined EIGEN_INTERNAL_DEBUGGING || defined EIGEN_BDCSVD_SANITY_CHECKS
RealScalar fRight = secularEq(rightShifted, col0, diag, perm, diagShifted, shift); RealScalar fRight = secularEq(rightShifted, col0, diag, perm, diagShifted, shift);
#endif #endif
#ifdef EIGEN_BDCSVD_SANITY_CHECKS #ifdef EIGEN_BDCSVD_SANITY_CHECKS
if(!(std::isfinite)(fLeft)) if(!(numext::isfinite)(fLeft))
std::cout << "f(" << leftShifted << ") =" << fLeft << " ; " << left << " " << shift << " " << right << "\n"; std::cout << "f(" << leftShifted << ") =" << fLeft << " ; " << left << " " << shift << " " << right << "\n";
assert((std::isfinite)(fLeft)); assert((numext::isfinite)(fLeft));
if(!(std::isfinite)(fRight)) if(!(numext::isfinite)(fRight))
std::cout << "f(" << rightShifted << ") =" << fRight << " ; " << left << " " << shift << " " << right << "\n"; std::cout << "f(" << rightShifted << ") =" << fRight << " ; " << left << " " << shift << " " << right << "\n";
// assert((std::isfinite)(fRight)); // assert((numext::isfinite)(fRight));
#endif #endif
#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
@ -885,29 +901,43 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
std::cout << "f(leftShifted) using leftShifted=" << leftShifted << " ; diagShifted(1:10):" << diagShifted.head(10).transpose() << "\n ; " std::cout << "f(leftShifted) using leftShifted=" << leftShifted << " ; diagShifted(1:10):" << diagShifted.head(10).transpose() << "\n ; "
<< "left==shift=" << bool(left==shift) << " ; left-shift = " << (left-shift) << "\n"; << "left==shift=" << bool(left==shift) << " ; left-shift = " << (left-shift) << "\n";
std::cout << "k=" << k << ", " << fLeft << " * " << fRight << " == " << fLeft * fRight << " ; " std::cout << "k=" << k << ", " << fLeft << " * " << fRight << " == " << fLeft * fRight << " ; "
<< "[" << left << " .. " << right << "] -> [" << leftShifted << " " << rightShifted << "], shift=" << shift << " , f(right)=" << secularEq(0, col0, diag, perm, diagShifted, shift) << " == " << secularEq(right, col0, diag, perm, diag, 0) << "\n"; << "[" << left << " .. " << right << "] -> [" << leftShifted << " " << rightShifted << "], shift=" << shift
<< " , f(right)=" << secularEq(0, col0, diag, perm, diagShifted, shift)
<< " == " << secularEq(right, col0, diag, perm, diag, 0) << " == " << fRight << "\n";
} }
#endif #endif
eigen_internal_assert(fLeft * fRight < Literal(0)); eigen_internal_assert(fLeft * fRight < Literal(0));
while (rightShifted - leftShifted > Literal(2) * NumTraits<RealScalar>::epsilon() * numext::maxi<RealScalar>(abs(leftShifted), abs(rightShifted))) if(fLeft<Literal(0))
{ {
RealScalar midShifted = (leftShifted + rightShifted) / Literal(2); while (rightShifted - leftShifted > Literal(2) * NumTraits<RealScalar>::epsilon() * numext::maxi<RealScalar>(abs(leftShifted), abs(rightShifted)))
fMid = secularEq(midShifted, col0, diag, perm, diagShifted, shift); {
eigen_internal_assert((numext::isfinite)(fMid)); RealScalar midShifted = (leftShifted + rightShifted) / Literal(2);
fMid = secularEq(midShifted, col0, diag, perm, diagShifted, shift);
eigen_internal_assert((numext::isfinite)(fMid));
if (fLeft * fMid < Literal(0)) if (fLeft * fMid < Literal(0))
{ {
rightShifted = midShifted; rightShifted = midShifted;
} }
else else
{ {
leftShifted = midShifted; leftShifted = midShifted;
fLeft = fMid; fLeft = fMid;
}
} }
muCur = (leftShifted + rightShifted) / Literal(2);
}
else
{
// We have a problem as shifting on the left or right give either a positive or negative value
// at the middle of [left,right]...
// Instead fo abbording or entering an infinite loop,
// let's just use the middle as the estimated zero-crossing:
muCur = (right - left) * RealScalar(0.5);
if(shift == right)
muCur = -muCur;
} }
muCur = (leftShifted + rightShifted) / Literal(2);
} }
singVals[k] = shift + muCur; singVals[k] = shift + muCur;
@ -994,7 +1024,7 @@ void BDCSVD<MatrixType>::perturbCol0
assert(prod>=0); assert(prod>=0);
#endif #endif
#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
if(i!=k && std::abs(((singVals(j)+dk)*(mus(j)+(shifts(j)-dk)))/((diag(i)+dk)*(diag(i)-dk)) - 1) > 0.9 ) if(i!=k && numext::abs(((singVals(j)+dk)*(mus(j)+(shifts(j)-dk)))/((diag(i)+dk)*(diag(i)-dk)) - 1) > 0.9 )
std::cout << " " << ((singVals(j)+dk)*(mus(j)+(shifts(j)-dk)))/((diag(i)+dk)*(diag(i)-dk)) << " == (" << (singVals(j)+dk) << " * " << (mus(j)+(shifts(j)-dk)) std::cout << " " << ((singVals(j)+dk)*(mus(j)+(shifts(j)-dk)))/((diag(i)+dk)*(diag(i)-dk)) << " == (" << (singVals(j)+dk) << " * " << (mus(j)+(shifts(j)-dk))
<< ") / (" << (diag(i)+dk) << " * " << (diag(i)-dk) << ")\n"; << ") / (" << (diag(i)+dk) << " * " << (diag(i)-dk) << ")\n";
#endif #endif
@ -1005,7 +1035,7 @@ void BDCSVD<MatrixType>::perturbCol0
#endif #endif
RealScalar tmp = sqrt(prod); RealScalar tmp = sqrt(prod);
#ifdef EIGEN_BDCSVD_SANITY_CHECKS #ifdef EIGEN_BDCSVD_SANITY_CHECKS
assert((std::isfinite)(tmp)); assert((numext::isfinite)(tmp));
#endif #endif
zhat(k) = col0(k) > Literal(0) ? RealScalar(tmp) : RealScalar(-tmp); zhat(k) = col0(k) > Literal(0) ? RealScalar(tmp) : RealScalar(-tmp);
} }