mirror of
https://gitlab.com/libeigen/eigen.git
synced 2025-08-11 03:09:01 +08:00
Improve debugging tests and output in BDCSVD
This commit is contained in:
Gael Guennebaud
parent
e8468ea91b
commit
e9d2888e74
@ -22,6 +22,11 @@
|
||||
// #define EIGEN_BDCSVD_DEBUG_VERBOSE
|
||||
// #define EIGEN_BDCSVD_SANITY_CHECKS
|
||||
|
||||
#ifdef EIGEN_BDCSVD_SANITY_CHECKS
|
||||
#undef eigen_internal_assert
|
||||
#define eigen_internal_assert(X) assert(X);
|
||||
#endif
|
||||
|
||||
namespace Eigen {
|
||||
|
||||
#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
|
||||
@ -591,7 +596,7 @@ void BDCSVD<MatrixType>::computeSVDofM(Index firstCol, Index n, MatrixXr& U, Vec
|
||||
// but others are interleaved and we must ignore them at this stage.
|
||||
// To this end, let's compute a permutation skipping them:
|
||||
Index actual_n = n;
|
||||
while(actual_n>1 && diag(actual_n-1)==Literal(0)) --actual_n;
|
||||
while(actual_n>1 && diag(actual_n-1)==Literal(0)) {--actual_n; eigen_internal_assert(col0(actual_n)==Literal(0)); }
|
||||
Index m = 0; // size of the deflated problem
|
||||
for(Index k=0;k<actual_n;++k)
|
||||
if(abs(col0(k))>considerZero)
|
||||
@ -618,13 +623,11 @@ void BDCSVD<MatrixType>::computeSVDofM(Index firstCol, Index n, MatrixXr& U, Vec
|
||||
std::cout << " shift: " << shifts.transpose() << "\n";
|
||||
|
||||
{
|
||||
Index actual_n = 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 << " check1 (expect0) : " << ((singVals.array()-(shifts+mus)) / singVals.array()).head(actual_n).transpose() << "\n\n";
|
||||
assert((((singVals.array()-(shifts+mus)) / singVals.array()).head(actual_n) >= 0).all());
|
||||
std::cout << " check2 (>0) : " << ((singVals.array()-diag) / singVals.array()).head(actual_n).transpose() << "\n\n";
|
||||
std::cout << " check3 (>0) : " << ((diag.segment(1,actual_n-1)-singVals.head(actual_n-1).array()) / singVals.head(actual_n-1).array()).transpose() << "\n\n\n";
|
||||
std::cout << " check4 (>0) : " << ((singVals.segment(1,actual_n-1)-singVals.head(actual_n-1))).transpose() << "\n\n\n";
|
||||
assert((((singVals.array()-diag) / singVals.array()).head(actual_n) >= 0).all());
|
||||
}
|
||||
#endif
|
||||
|
||||
@ -652,13 +655,13 @@ void BDCSVD<MatrixType>::computeSVDofM(Index firstCol, Index n, MatrixXr& U, Vec
|
||||
#endif
|
||||
|
||||
#ifdef EIGEN_BDCSVD_SANITY_CHECKS
|
||||
assert(U.allFinite());
|
||||
assert(V.allFinite());
|
||||
assert((U.transpose() * U - MatrixXr(MatrixXr::Identity(U.cols(),U.cols()))).norm() < 1e-14 * n);
|
||||
assert((V.transpose() * V - MatrixXr(MatrixXr::Identity(V.cols(),V.cols()))).norm() < 1e-14 * n);
|
||||
assert(m_naiveU.allFinite());
|
||||
assert(m_naiveV.allFinite());
|
||||
assert(m_computed.allFinite());
|
||||
assert(U.allFinite());
|
||||
assert(V.allFinite());
|
||||
// assert((U.transpose() * U - MatrixXr(MatrixXr::Identity(U.cols(),U.cols()))).norm() < 100*NumTraits<RealScalar>::epsilon() * n);
|
||||
// assert((V.transpose() * V - MatrixXr(MatrixXr::Identity(V.cols(),V.cols()))).norm() < 100*NumTraits<RealScalar>::epsilon() * n);
|
||||
#endif
|
||||
|
||||
// Because of deflation, the singular values might not be completely sorted.
|
||||
@ -673,6 +676,15 @@ void BDCSVD<MatrixType>::computeSVDofM(Index firstCol, Index n, MatrixXr& U, Vec
|
||||
if(m_compV) V.col(i).swap(V.col(i+1));
|
||||
}
|
||||
}
|
||||
|
||||
#ifdef EIGEN_BDCSVD_SANITY_CHECKS
|
||||
{
|
||||
bool singular_values_sorted = (((singVals.segment(1,actual_n-1)-singVals.head(actual_n-1))).array() >= 0).all();
|
||||
if(!singular_values_sorted)
|
||||
std::cout << "Singular values are not sorted: " << singVals.segment(1,actual_n).transpose() << "\n";
|
||||
assert(singular_values_sorted);
|
||||
}
|
||||
#endif
|
||||
|
||||
// Reverse order so that singular values in increased order
|
||||
// Because of deflation, the zeros singular-values are already at the end
|
||||
@ -749,19 +761,22 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
|
||||
RealScalar mid = left + (right-left) / Literal(2);
|
||||
RealScalar fMid = secularEq(mid, col0, diag, perm, diag, Literal(0));
|
||||
#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
|
||||
std::cout << right-left << "\n";
|
||||
std::cout << "fMid = " << fMid << " " << secularEq(mid-left, col0, diag, perm, diag-left, left) << " " << secularEq(mid-right, col0, diag, perm, diag-right, right) << "\n";
|
||||
std::cout << " = " << secularEq(0.1*(left+right), col0, diag, perm, diag, 0)
|
||||
<< " " << secularEq(0.2*(left+right), col0, diag, perm, diag, 0)
|
||||
<< " " << secularEq(0.3*(left+right), col0, diag, perm, diag, 0)
|
||||
<< " " << secularEq(0.4*(left+right), col0, diag, perm, diag, 0)
|
||||
<< " " << secularEq(0.49*(left+right), col0, diag, perm, diag, 0)
|
||||
<< " " << secularEq(0.5*(left+right), col0, diag, perm, diag, 0)
|
||||
<< " " << secularEq(0.51*(left+right), col0, diag, perm, diag, 0)
|
||||
<< " " << secularEq(0.6*(left+right), col0, diag, perm, diag, 0)
|
||||
<< " " << secularEq(0.7*(left+right), col0, diag, perm, diag, 0)
|
||||
<< " " << secularEq(0.8*(left+right), col0, diag, perm, diag, 0)
|
||||
<< " " << secularEq(0.9*(left+right), col0, diag, perm, diag, 0) << "\n";
|
||||
std::cout << "right-left = " << right-left << "\n";
|
||||
// std::cout << "fMid = " << fMid << " " << secularEq(mid-left, col0, diag, perm, ArrayXr(diag-left), left)
|
||||
// << " " << secularEq(mid-right, col0, diag, perm, ArrayXr(diag-right), right) << "\n";
|
||||
std::cout << " = " << secularEq(left+RealScalar(0.000001)*(right-left), col0, diag, perm, diag, 0)
|
||||
<< " " << secularEq(left+RealScalar(0.1) *(right-left), col0, diag, perm, diag, 0)
|
||||
<< " " << secularEq(left+RealScalar(0.2) *(right-left), col0, diag, perm, diag, 0)
|
||||
<< " " << secularEq(left+RealScalar(0.3) *(right-left), col0, diag, perm, diag, 0)
|
||||
<< " " << secularEq(left+RealScalar(0.4) *(right-left), col0, diag, perm, diag, 0)
|
||||
<< " " << secularEq(left+RealScalar(0.49) *(right-left), col0, diag, perm, diag, 0)
|
||||
<< " " << secularEq(left+RealScalar(0.5) *(right-left), col0, diag, perm, diag, 0)
|
||||
<< " " << secularEq(left+RealScalar(0.51) *(right-left), col0, diag, perm, diag, 0)
|
||||
<< " " << secularEq(left+RealScalar(0.6) *(right-left), col0, diag, perm, diag, 0)
|
||||
<< " " << secularEq(left+RealScalar(0.7) *(right-left), col0, diag, perm, diag, 0)
|
||||
<< " " << secularEq(left+RealScalar(0.8) *(right-left), col0, diag, perm, diag, 0)
|
||||
<< " " << secularEq(left+RealScalar(0.9) *(right-left), col0, diag, perm, diag, 0)
|
||||
<< " " << secularEq(left+RealScalar(0.999999)*(right-left), col0, diag, perm, diag, 0) << "\n";
|
||||
#endif
|
||||
RealScalar shift = (k == actual_n-1 || fMid > Literal(0)) ? left : right;
|
||||
|
||||
@ -804,13 +819,16 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
|
||||
// And find mu such that f(mu)==0:
|
||||
RealScalar muZero = -a/b;
|
||||
RealScalar fZero = secularEq(muZero, col0, diag, perm, diagShifted, shift);
|
||||
|
||||
#ifdef EIGEN_BDCSVD_SANITY_CHECKS
|
||||
assert((std::isfinite)(fZero));
|
||||
#endif
|
||||
|
||||
muPrev = muCur;
|
||||
fPrev = fCur;
|
||||
muCur = muZero;
|
||||
fCur = fZero;
|
||||
|
||||
|
||||
if (shift == left && (muCur < Literal(0) || muCur > right - left)) useBisection = true;
|
||||
if (shift == right && (muCur < -(right - left) || muCur > Literal(0))) useBisection = true;
|
||||
if (abs(fCur)>abs(fPrev)) useBisection = true;
|
||||
@ -843,18 +861,30 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
|
||||
else
|
||||
rightShifted = -(std::numeric_limits<RealScalar>::min)();
|
||||
}
|
||||
|
||||
|
||||
RealScalar fLeft = secularEq(leftShifted, col0, diag, perm, diagShifted, shift);
|
||||
|
||||
#if defined EIGEN_INTERNAL_DEBUGGING || defined EIGEN_BDCSVD_DEBUG_VERBOSE
|
||||
#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))
|
||||
std::cout << "f(" << leftShifted << ") =" << fLeft << " ; " << left << " " << shift << " " << right << "\n";
|
||||
assert((std::isfinite)(fLeft));
|
||||
|
||||
if(!(std::isfinite)(fRight))
|
||||
std::cout << "f(" << rightShifted << ") =" << fRight << " ; " << left << " " << shift << " " << right << "\n";
|
||||
// assert((std::isfinite)(fRight));
|
||||
#endif
|
||||
|
||||
#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
|
||||
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 << "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";
|
||||
}
|
||||
#endif
|
||||
eigen_internal_assert(fLeft * fRight < Literal(0));
|
||||
@ -863,6 +893,8 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
|
||||
{
|
||||
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;
|
||||
@ -881,6 +913,15 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
|
||||
shifts[k] = shift;
|
||||
mus[k] = muCur;
|
||||
|
||||
#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
|
||||
if(k+1<n)
|
||||
std::cout << "found " << singVals[k] << " == " << shift << " + " << muCur << " from " << diag(k) << " .. " << diag(k+1) << "\n";
|
||||
#endif
|
||||
#ifdef EIGEN_BDCSVD_SANITY_CHECKS
|
||||
assert(k==0 || singVals[k]>=singVals[k-1]);
|
||||
assert(singVals[k]>=diag(k));
|
||||
#endif
|
||||
|
||||
// perturb singular value slightly if it equals diagonal entry to avoid division by zero later
|
||||
// (deflation is supposed to avoid this from happening)
|
||||
// - this does no seem to be necessary anymore -
|
||||
@ -915,14 +956,42 @@ void BDCSVD<MatrixType>::perturbCol0
|
||||
// see equation (3.6)
|
||||
RealScalar dk = diag(k);
|
||||
RealScalar prod = (singVals(last) + dk) * (mus(last) + (shifts(last) - dk));
|
||||
#ifdef EIGEN_BDCSVD_SANITY_CHECKS
|
||||
if(prod<0) {
|
||||
std::cout << "k = " << k << " ; z(k)=" << col0(k) << ", diag(k)=" << dk << "\n";
|
||||
std::cout << "prod = " << "(" << singVals(last) << " + " << dk << ") * (" << mus(last) << " + (" << shifts(last) << " - " << dk << "))" << "\n";
|
||||
std::cout << " = " << singVals(last) + dk << " * " << mus(last) + (shifts(last) - dk) << "\n";
|
||||
}
|
||||
assert(prod>=0);
|
||||
#endif
|
||||
|
||||
for(Index l = 0; l<m; ++l)
|
||||
{
|
||||
Index i = perm(l);
|
||||
if(i!=k)
|
||||
{
|
||||
#ifdef EIGEN_BDCSVD_SANITY_CHECKS
|
||||
if(i>=k && (l==0 || l-1>=m))
|
||||
{
|
||||
std::cout << "Error in perturbCol0\n";
|
||||
std::cout << " " << k << "/" << n << " " << l << "/" << m << " " << i << "/" << n << " ; " << col0(k) << " " << diag(k) << " " << "\n";
|
||||
std::cout << " " <<diag(i) << "\n";
|
||||
Index j = (i<k /*|| l==0*/) ? i : perm(l-1);
|
||||
std::cout << " " << "j=" << j << "\n";
|
||||
}
|
||||
#endif
|
||||
Index j = i<k ? i : perm(l-1);
|
||||
#ifdef EIGEN_BDCSVD_SANITY_CHECKS
|
||||
if(!(dk!=Literal(0) || diag(i)!=Literal(0)))
|
||||
{
|
||||
std::cout << "k=" << k << ", i=" << i << ", l=" << l << ", perm.size()=" << perm.size() << "\n";
|
||||
}
|
||||
assert(dk!=Literal(0) || diag(i)!=Literal(0));
|
||||
#endif
|
||||
prod *= ((singVals(j)+dk) / ((diag(i)+dk))) * ((mus(j)+(shifts(j)-dk)) / ((diag(i)-dk)));
|
||||
#ifdef EIGEN_BDCSVD_SANITY_CHECKS
|
||||
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 )
|
||||
std::cout << " " << ((singVals(j)+dk)*(mus(j)+(shifts(j)-dk)))/((diag(i)+dk)*(diag(i)-dk)) << " == (" << (singVals(j)+dk) << " * " << (mus(j)+(shifts(j)-dk))
|
||||
@ -934,6 +1003,9 @@ void BDCSVD<MatrixType>::perturbCol0
|
||||
std::cout << "zhat(" << k << ") = sqrt( " << prod << ") ; " << (singVals(last) + dk) << " * " << mus(last) + shifts(last) << " - " << dk << "\n";
|
||||
#endif
|
||||
RealScalar tmp = sqrt(prod);
|
||||
#ifdef EIGEN_BDCSVD_SANITY_CHECKS
|
||||
assert((std::isfinite)(tmp));
|
||||
#endif
|
||||
zhat(k) = col0(k) > Literal(0) ? tmp : -tmp;
|
||||
}
|
||||
}
|
||||
@ -1043,7 +1115,7 @@ void BDCSVD<MatrixType>::deflation44(Index firstColu , Index firstColm, Index fi
|
||||
}
|
||||
c/=r;
|
||||
s/=r;
|
||||
m_computed(firstColm + i, firstColm) = r;
|
||||
m_computed(firstColm + i, firstColm) = r;
|
||||
m_computed(firstColm + j, firstColm + j) = m_computed(firstColm + i, firstColm + i);
|
||||
m_computed(firstColm + j, firstColm) = Literal(0);
|
||||
|
||||
@ -1117,6 +1189,7 @@ void BDCSVD<MatrixType>::deflation(Index firstCol, Index lastCol, Index k, Index
|
||||
#endif
|
||||
#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
|
||||
std::cout << "to be sorted: " << diag.transpose() << "\n\n";
|
||||
std::cout << " : " << col0.transpose() << "\n\n";
|
||||
#endif
|
||||
{
|
||||
// Check for total deflation
|
||||
@ -1207,7 +1280,7 @@ void BDCSVD<MatrixType>::deflation(Index firstCol, Index lastCol, Index k, Index
|
||||
if( (diag(i) - diag(i-1)) < NumTraits<RealScalar>::epsilon()*maxDiag )
|
||||
{
|
||||
#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
|
||||
std::cout << "deflation 4.4 with i = " << i << " because " << (diag(i) - diag(i-1)) << " < " << NumTraits<RealScalar>::epsilon()*diag(i) << "\n";
|
||||
std::cout << "deflation 4.4 with i = " << i << " because " << diag(i) << " - " << diag(i-1) << " == " << (diag(i) - diag(i-1)) << " < " << NumTraits<RealScalar>::epsilon()*/*diag(i)*/maxDiag << "\n";
|
||||
#endif
|
||||
eigen_internal_assert(abs(diag(i) - diag(i-1))<epsilon_coarse && " diagonal entries are not properly sorted");
|
||||
deflation44(firstCol, firstCol + shift, firstRowW, firstColW, i-1, i, length);
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user