diff options
author | Gael Guennebaud <g.gael@free.fr> | 2019-03-27 20:16:58 +0100 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2019-03-27 20:16:58 +0100 |
commit | 45e65fbb7791e453f88f959111cff45e0fb7dd6f (patch) | |
tree | 995c809517d9cb4493ed83a8f08ed46fb6af65e3 /Eigen/src/SVD | |
parent | 8de66719f91dd2b3ed621e174715c46c97b63ce4 (diff) |
bug #1695: fix a numerical robustness issue. Computing the secular equation at the middle range without a shift might give a wrong sign.
Diffstat (limited to 'Eigen/src/SVD')
-rw-r--r-- | Eigen/src/SVD/BDCSVD.h | 78 |
1 files changed, 54 insertions, 24 deletions
diff --git a/Eigen/src/SVD/BDCSVD.h b/Eigen/src/SVD/BDCSVD.h index e3fddacbc..bcec45f58 100644 --- a/Eigen/src/SVD/BDCSVD.h +++ b/Eigen/src/SVD/BDCSVD.h @@ -784,6 +784,21 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d // measure everything relative to shift Map<ArrayXr> diagShifted(m_workspace.data()+4*n, n); 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 RealScalar muPrev, muCur; @@ -822,7 +837,7 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d RealScalar fZero = secularEq(muZero, col0, diag, perm, diagShifted, shift); #ifdef EIGEN_BDCSVD_SANITY_CHECKS - assert((std::isfinite)(fZero)); + assert((numext::isfinite)(fZero)); #endif muPrev = muCur; @@ -864,50 +879,65 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d } 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 RealScalar fRight = secularEq(rightShifted, col0, diag, perm, diagShifted, shift); #endif #ifdef EIGEN_BDCSVD_SANITY_CHECKS - if(!(std::isfinite)(fLeft)) + if(!(numext::isfinite)(fLeft)) 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"; -// assert((std::isfinite)(fRight)); + // assert((numext::isfinite)(fRight)); #endif - + #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE if(!(fLeft * fRight<0)) { 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"; 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 eigen_internal_assert(fLeft * fRight < Literal(0)); - - while (rightShifted - leftShifted > Literal(2) * NumTraits<RealScalar>::epsilon() * numext::maxi<RealScalar>(abs(leftShifted), abs(rightShifted))) - { - 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)) - { - rightShifted = midShifted; - } - else + if(fLeft<Literal(0)) + { + while (rightShifted - leftShifted > Literal(2) * NumTraits<RealScalar>::epsilon() * numext::maxi<RealScalar>(abs(leftShifted), abs(rightShifted))) { - leftShifted = midShifted; - fLeft = 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)) + { + rightShifted = midShifted; + } + else + { + leftShifted = midShifted; + 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; @@ -994,7 +1024,7 @@ void BDCSVD<MatrixType>::perturbCol0 assert(prod>=0); #endif #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)) << ") / (" << (diag(i)+dk) << " * " << (diag(i)-dk) << ")\n"; #endif @@ -1005,7 +1035,7 @@ void BDCSVD<MatrixType>::perturbCol0 #endif RealScalar tmp = sqrt(prod); #ifdef EIGEN_BDCSVD_SANITY_CHECKS - assert((std::isfinite)(tmp)); + assert((numext::isfinite)(tmp)); #endif zhat(k) = col0(k) > Literal(0) ? RealScalar(tmp) : RealScalar(-tmp); } |