From 2cc41dbe8355eee2e646c5bffa2d4b6cfdad4029 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Wed, 15 Oct 2014 15:21:12 +0200 Subject: [PATCH] D&C SVD: fix some numerical issues by truly skipping deflated singular values when computing them --- unsupported/Eigen/src/BDCSVD/BDCSVD.h | 143 +++++++++++--------------- 1 file changed, 59 insertions(+), 84 deletions(-) diff --git a/unsupported/Eigen/src/BDCSVD/BDCSVD.h b/unsupported/Eigen/src/BDCSVD/BDCSVD.h index d8d75624d..e8e795589 100644 --- a/unsupported/Eigen/src/BDCSVD/BDCSVD.h +++ b/unsupported/Eigen/src/BDCSVD/BDCSVD.h @@ -84,6 +84,7 @@ public: typedef Matrix MatrixXr; typedef Matrix VectorType; typedef Array ArrayXr; + typedef Array ArrayXi; /** \brief Default Constructor. * @@ -159,19 +160,16 @@ private: void allocate(Index rows, Index cols, unsigned int computationOptions); void divide(Index firstCol, Index lastCol, Index firstRowW, Index firstColW, Index shift); void computeSVDofM(Index firstCol, Index n, MatrixXr& U, VectorType& singVals, MatrixXr& V); - void computeSingVals(const ArrayXr& col0, const ArrayXr& diag, VectorType& singVals, - ArrayXr& shifts, ArrayXr& mus); - void perturbCol0(const ArrayXr& col0, const ArrayXr& diag, const VectorType& singVals, - const ArrayXr& shifts, const ArrayXr& mus, ArrayXr& zhat); - void computeSingVecs(const ArrayXr& zhat, const ArrayXr& diag, const VectorType& singVals, - const ArrayXr& shifts, const ArrayXr& mus, MatrixXr& U, MatrixXr& V); + void computeSingVals(const ArrayXr& col0, const ArrayXr& diag, const ArrayXi& perm, VectorType& singVals, ArrayXr& shifts, ArrayXr& mus); + void perturbCol0(const ArrayXr& col0, const ArrayXr& diag, const ArrayXi& perm, const VectorType& singVals, const ArrayXr& shifts, const ArrayXr& mus, ArrayXr& zhat); + void computeSingVecs(const ArrayXr& zhat, const ArrayXr& diag, const ArrayXi& perm, const VectorType& singVals, const ArrayXr& shifts, const ArrayXr& mus, MatrixXr& U, MatrixXr& V); void deflation43(Index firstCol, Index shift, Index i, Index size); void deflation44(Index firstColu , Index firstColm, Index firstRowW, Index firstColW, Index i, Index j, Index size); void deflation(Index firstCol, Index lastCol, Index k, Index firstRowW, Index firstColW, Index shift); template void copyUV(const HouseholderU &householderU, const HouseholderV &householderV, const NaiveU &naiveU, const NaiveV &naivev); static void structured_update(Block A, const MatrixXr &B, Index n1); - static RealScalar secularEq(RealScalar x, const ArrayXr& col0, const ArrayXr& diag, const ArrayXr& diagShifted, RealScalar shift, Index n); + static RealScalar secularEq(RealScalar x, const ArrayXr& col0, const ArrayXr& diag, const ArrayXi &perm, const ArrayXr& diagShifted, RealScalar shift); protected: MatrixXr m_naiveU, m_naiveV; @@ -543,13 +541,24 @@ void BDCSVD::computeSVDofM(Index firstCol, Index n, MatrixXr& U, Vec ArrayXr diag = m_computed.block(firstCol, firstCol, n, n).diagonal(); diag(0) = 0; - // compute singular values and vectors (in decreasing order) + // Allocate space for singular values and vectors singVals.resize(n); U.resize(n+1, n+1); if (m_compV) V.resize(n, n); if (col0.hasNaN() || diag.hasNaN()) { std::cout << "\n\nHAS NAN\n\n"; return; } - + + // Many singular values might have been deflated, the zero ones have been moved to the end, + // 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)==0) --actual_n; + Index m = 0; // size of the deflated problem + ArrayXi perm(actual_n); + for(Index k=0;k::computeSVDofM(Index firstCol, Index n, MatrixXr& U, Vec #endif // Compute singVals, shifts, and mus - computeSingVals(col0, diag, singVals, shifts, mus); + computeSingVals(col0, diag, perm, singVals, shifts, mus); #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE std::cout << " j: " << (m_computed.block(firstCol, firstCol, n, n)).jacobiSvd().singularValues().transpose().reverse() << "\n\n"; @@ -586,7 +595,7 @@ void BDCSVD::computeSVDofM(Index firstCol, Index n, MatrixXr& U, Vec #endif // Compute zhat - perturbCol0(col0, diag, singVals, shifts, mus, zhat); + perturbCol0(col0, diag, perm, singVals, shifts, mus, zhat); #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE std::cout << " zhat: " << zhat.transpose() << "\n"; #endif @@ -595,7 +604,7 @@ void BDCSVD::computeSVDofM(Index firstCol, Index n, MatrixXr& U, Vec assert(zhat.allFinite()); #endif - computeSingVecs(zhat, diag, singVals, shifts, mus, U, V); + computeSingVecs(zhat, diag, perm, singVals, shifts, mus, U, V); #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE std::cout << "U^T U: " << (U.transpose() * U - MatrixXr(MatrixXr::Identity(U.cols(),U.cols()))).norm() << "\n"; @@ -610,9 +619,6 @@ void BDCSVD::computeSVDofM(Index firstCol, Index n, MatrixXr& U, Vec assert(m_computed.allFinite()); #endif - Index actual_n = n; - while(actual_n>1 && singVals(actual_n-1)==0) --actual_n; - // Because of deflation, the singular values might not be completely sorted. // Fortunately, reordering them is a O(n) problem for(Index i=0; i::computeSVDofM(Index firstCol, Index n, MatrixXr& U, Vec } template -typename BDCSVD::RealScalar BDCSVD::secularEq(RealScalar mu, const ArrayXr& col0, const ArrayXr& diag, const ArrayXr& diagShifted, RealScalar shift, Index n) +typename BDCSVD::RealScalar BDCSVD::secularEq(RealScalar mu, const ArrayXr& col0, const ArrayXr& diag, const ArrayXi &perm, const ArrayXr& diagShifted, RealScalar shift) { - return 1 + (col0.square() / ((diagShifted - mu) * (diag + shift + mu))).head(n).sum(); + Index m = perm.size(); + RealScalar res = 1; + for(Index i=0; i -void BDCSVD::computeSingVals(const ArrayXr& col0, const ArrayXr& diag, +void BDCSVD::computeSingVals(const ArrayXr& col0, const ArrayXr& diag, const ArrayXi &perm, VectorType& singVals, ArrayXr& shifts, ArrayXr& mus) { using std::abs; @@ -657,16 +670,6 @@ void BDCSVD::computeSingVals(const ArrayXr& col0, const ArrayXr& dia Index n = col0.size(); Index actual_n = n; while(actual_n>1 && col0(actual_n-1)==0) --actual_n; -// Index m = 0; -// Array perm(actual_n); -// { -// for(Index k=0;k::computeSingVals(const ArrayXr& col0, const ArrayXr& dia Index l = k+1; while(col0(l)==0) { ++l; eigen_internal_assert(l 0) ? left : right; @@ -727,8 +733,8 @@ void BDCSVD::computeSingVals(const ArrayXr& col0, const ArrayXr& dia muCur = -(right - left) * 0.5; } - RealScalar fPrev = secularEq(muPrev, col0, diag, diagShifted, shift, actual_n); - RealScalar fCur = secularEq(muCur, col0, diag, diagShifted, shift, actual_n); + RealScalar fPrev = secularEq(muPrev, col0, diag, perm, diagShifted, shift); + RealScalar fCur = secularEq(muCur, col0, diag, perm, diagShifted, shift); if (abs(fPrev) < abs(fCur)) { swap(fPrev, fCur); @@ -747,7 +753,7 @@ void BDCSVD::computeSingVals(const ArrayXr& col0, const ArrayXr& dia RealScalar b = fCur - a / muCur; // And find mu such that f(mu)==0: RealScalar muZero = -a/b; - RealScalar fZero = secularEq(muZero, col0, diag, diagShifted, shift, actual_n); + RealScalar fZero = secularEq(muZero, col0, diag, perm, diagShifted, shift); muPrev = muCur; fPrev = fCur; @@ -780,8 +786,8 @@ void BDCSVD::computeSingVals(const ArrayXr& col0, const ArrayXr& dia rightShifted = -RealScalar(1)/NumTraits::highest(); } - RealScalar fLeft = secularEq(leftShifted, col0, diag, diagShifted, shift, actual_n); - RealScalar fRight = secularEq(rightShifted, col0, diag, diagShifted, shift, actual_n); + RealScalar fLeft = secularEq(leftShifted, col0, diag, perm, diagShifted, shift); + RealScalar fRight = secularEq(rightShifted, col0, diag, perm, diagShifted, shift); #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE if(!(fLeft * fRight<0)) @@ -792,7 +798,7 @@ void BDCSVD::computeSingVals(const ArrayXr& col0, const ArrayXr& dia while (rightShifted - leftShifted > 2 * NumTraits::epsilon() * (max)(abs(leftShifted), abs(rightShifted))) { RealScalar midShifted = (leftShifted + rightShifted) / 2; - RealScalar fMid = secularEq(midShifted, col0, diag, diagShifted, shift, actual_n); + RealScalar fMid = secularEq(midShifted, col0, diag, perm, diagShifted, shift); if (fLeft * fMid < 0) { rightShifted = midShifted; @@ -824,28 +830,13 @@ void BDCSVD::computeSingVals(const ArrayXr& col0, const ArrayXr& dia // zhat is perturbation of col0 for which singular vectors can be computed stably (see Section 3.1) template void BDCSVD::perturbCol0 - (const ArrayXr& col0, const ArrayXr& diag, const VectorType& singVals, + (const ArrayXr& col0, const ArrayXr& diag, const ArrayXi &perm, const VectorType& singVals, const ArrayXr& shifts, const ArrayXr& mus, ArrayXr& zhat) { using std::sqrt; Index n = col0.size(); - - // Ignore trailing zeros: - Index actual_n = n; - while(actual_n>1 && col0(actual_n-1)==0) --actual_n; - // Deflated non-zero singular values are kept in-place, - // we thus compute an indirection array to properly ignore all deflated entries. - // TODO compute it once! - Index m = 0; // size of the deflated problem - Array perm(actual_n); - { - for(Index k=0;k::perturbCol0 { // see equation (3.6) RealScalar dk = diag(k); - RealScalar prod = (singVals(actual_n-1) + dk) * (mus(actual_n-1) + (shifts(actual_n-1) - dk)); + RealScalar prod = (singVals(last) + dk) * (mus(last) + (shifts(last) - dk)); for(Index l = 0; l::perturbCol0 } } #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE - std::cout << "zhat(" << k << ") = sqrt( " << prod << ") ; " << (singVals(actual_n-1) + dk) << " * " << mus(actual_n-1) + shifts(actual_n-1) << " - " << dk << "\n"; + std::cout << "zhat(" << k << ") = sqrt( " << prod << ") ; " << (singVals(last) + dk) << " * " << mus(last) + shifts(last) << " - " << dk << "\n"; #endif RealScalar tmp = sqrt(prod); zhat(k) = col0(k) > 0 ? tmp : -tmp; @@ -884,26 +875,11 @@ void BDCSVD::perturbCol0 // compute singular vectors template void BDCSVD::computeSingVecs - (const ArrayXr& zhat, const ArrayXr& diag, const VectorType& singVals, + (const ArrayXr& zhat, const ArrayXr& diag, const ArrayXi &perm, const VectorType& singVals, const ArrayXr& shifts, const ArrayXr& mus, MatrixXr& U, MatrixXr& V) { Index n = zhat.size(); - - // Deflated non-zero singular values are kept in-place, - // we thus compute an indirection array to properly ignore all deflated entries. - // TODO compute it once! - Index actual_n = n; - while(actual_n>1 && zhat(actual_n-1)==0) --actual_n; - Index m = 0; - Array perm(actual_n); - { - for(Index k=0;k::deflation44(Index firstColu , Index firstColm, Index fi c/=r; s/=r; m_computed(firstColm + i, firstColm) = r; - m_computed(firstColm + i, firstColm + i) = m_computed(firstColm + j, firstColm + j); + m_computed(firstColm + j, firstColm + j) = m_computed(firstColm + i, firstColm + i); m_computed(firstColm + j, firstColm) = 0; JacobiRotation J(c,-s); @@ -1117,7 +1093,6 @@ void BDCSVD::deflation(Index firstCol, Index lastCol, Index k, Index } } } -// std::cout << "perm: " << Matrix::Map(permutation, length) << "\n\n"; // Current index of each col, and current column of each index Index *realInd = new Index[length]; // FIXME avoid repeated dynamic memory allocation @@ -1165,7 +1140,7 @@ void BDCSVD::deflation(Index firstCol, Index lastCol, Index k, Index Index i = length-1; while(i>0 && (diag(i)==0 || col0(i)==0)) --i; for(; i>1;--i) - if( (diag(i) - diag(i-1)) < NumTraits::epsilon()*diag(i) ) + if( (diag(i) - diag(i-1)) < NumTraits::epsilon()*maxDiag ) { #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE std::cout << "deflation 4.4 with i = " << i << " because " << (diag(i) - diag(i-1)) << " < " << NumTraits::epsilon()*diag(i) << "\n";