merge
diff --git a/test/svd_common.h b/test/svd_common.h index e902d23..347ea80 100644 --- a/test/svd_common.h +++ b/test/svd_common.h
@@ -38,7 +38,6 @@ sigma.diagonal() = svd.singularValues().template cast<Scalar>(); MatrixUType u = svd.matrixU(); MatrixVType v = svd.matrixV(); - VERIFY_IS_APPROX(m, u * sigma * v.adjoint()); VERIFY_IS_UNITARY(u); VERIFY_IS_UNITARY(v); @@ -90,31 +89,31 @@ SolutionType x = svd.solve(rhs); + // evaluate normal equation which works also for least-squares solutions + if(internal::is_same<RealScalar,double>::value || svd.rank()==m.diagonal().size()) + { + // This test is not stable with single precision. + // This is probably because squaring m signicantly affects the precision. + VERIFY_IS_APPROX(m.adjoint()*(m*x),m.adjoint()*rhs); + } + RealScalar residual = (m*x-rhs).norm(); // Check that there is no significantly better solution in the neighborhood of x if(!test_isMuchSmallerThan(residual,rhs.norm())) { - // If the residual is very small, then we have an exact solution, so we are already good. - for(int k=0;k<x.rows();++k) + // ^^^ If the residual is very small, then we have an exact solution, so we are already good. + for(Index k=0;k<x.rows();++k) { SolutionType y(x); - y.row(k).array() += 2*NumTraits<RealScalar>::epsilon(); + y.row(k) = (1.+2*NumTraits<RealScalar>::epsilon())*x.row(k); RealScalar residual_y = (m*y-rhs).norm(); VERIFY( test_isApprox(residual_y,residual) || residual < residual_y ); - y.row(k) = x.row(k).array() - 2*NumTraits<RealScalar>::epsilon(); + y.row(k) = (1.-2*NumTraits<RealScalar>::epsilon())*x.row(k); residual_y = (m*y-rhs).norm(); VERIFY( test_isApprox(residual_y,residual) || residual < residual_y ); } } - - // evaluate normal equation which works also for least-squares solutions - if(internal::is_same<RealScalar,double>::value) - { - // This test is not stable with single precision. - // This is probably because squaring m signicantly affects the precision. - VERIFY_IS_APPROX(m.adjoint()*m*x,m.adjoint()*rhs); - } } // check minimal norm solutions, the inoput matrix m is only used to recover problem size @@ -234,11 +233,49 @@ Matrix<RealScalar,Dynamic,1> d = Matrix<RealScalar,Dynamic,1>::Random(diagSize); for(Index k=0; k<diagSize; ++k) d(k) = d(k)*std::pow(RealScalar(10),internal::random<RealScalar>(-s,s)); - m = Matrix<Scalar,Dynamic,Dynamic>::Random(m.rows(),diagSize) * d.asDiagonal() * Matrix<Scalar,Dynamic,Dynamic>::Random(diagSize,m.cols()); + + bool dup = internal::random<int>(0,10) < 3; + bool unit_uv = internal::random<int>(0,10) < (dup?7:3); // if we duplicate some diagonal entries, then increase the chance to preserve them using unitary U and V factors + + // duplicate some singular values + if(dup) + { + Index n = internal::random<Index>(0,d.size()-1); + for(Index i=0; i<n; ++i) + d(internal::random<Index>(0,d.size()-1)) = d(internal::random<Index>(0,d.size()-1)); + } + + Matrix<Scalar,Dynamic,Dynamic> U(m.rows(),diagSize); + Matrix<Scalar,Dynamic,Dynamic> VT(diagSize,m.cols()); + if(unit_uv) + { + // in very rare cases let's try with a pure diagonal matrix + if(internal::random<int>(0,10) < 1) + { + U.setIdentity(); + VT.setIdentity(); + } + else + { + createRandomPIMatrixOfRank(diagSize,U.rows(), U.cols(), U); + createRandomPIMatrixOfRank(diagSize,VT.rows(), VT.cols(), VT); + } + } + else + { + U.setRandom(); + VT.setRandom(); + } + + m = U * d.asDiagonal() * VT; + // cancel some coeffs - Index n = internal::random<Index>(0,m.size()-1); - for(Index i=0; i<n; ++i) - m(internal::random<Index>(0,m.rows()-1), internal::random<Index>(0,m.cols()-1)) = Scalar(0); + if(!(dup && unit_uv)) + { + Index n = internal::random<Index>(0,m.size()-1); + for(Index i=0; i<n; ++i) + m(internal::random<Index>(0,m.rows()-1), internal::random<Index>(0,m.cols()-1)) = Scalar(0); + } }
diff --git a/unsupported/Eigen/src/BDCSVD/BDCSVD.h b/unsupported/Eigen/src/BDCSVD/BDCSVD.h index d5e8140..e8e7955 100644 --- a/unsupported/Eigen/src/BDCSVD/BDCSVD.h +++ b/unsupported/Eigen/src/BDCSVD/BDCSVD.h
@@ -23,6 +23,10 @@ // #define EIGEN_BDCSVD_SANITY_CHECKS namespace Eigen { +#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE +IOFormat bdcsvdfmt(8, 0, ", ", "\n", " [", "]"); +#endif + template<typename _MatrixType> class BDCSVD; namespace internal { @@ -80,6 +84,7 @@ typedef Matrix<RealScalar, Dynamic, Dynamic> MatrixXr; typedef Matrix<RealScalar, Dynamic, 1> VectorType; typedef Array<RealScalar, Dynamic, 1> ArrayXr; + typedef Array<Index,1,Dynamic> ArrayXi; /** \brief Default Constructor. * @@ -155,19 +160,16 @@ 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<typename HouseholderU, typename HouseholderV, typename NaiveU, typename NaiveV> void copyUV(const HouseholderU &householderU, const HouseholderV &householderV, const NaiveU &naiveU, const NaiveV &naivev); static void structured_update(Block<MatrixXr,Dynamic,Dynamic> 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; @@ -234,6 +236,9 @@ template<typename MatrixType> BDCSVD<MatrixType>& BDCSVD<MatrixType>::compute(const MatrixType& matrix, unsigned int computationOptions) { +#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE + std::cout << "\n\n\n======================================================================================================================\n\n\n"; +#endif allocate(matrix.rows(), matrix.cols(), computationOptions); using std::abs; @@ -478,9 +483,23 @@ m_computed.col(firstCol + shift).segment(firstCol + shift + 1, k) = alphaK * l.transpose().real(); m_computed.col(firstCol + shift).segment(firstCol + shift + k + 1, n - k - 1) = betaK * f.transpose().real(); +#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE + ArrayXr tmp1 = (m_computed.block(firstCol+shift, firstCol+shift, n, n)).jacobiSvd().singularValues(); +#endif // Second part: try to deflate singular values in combined matrix deflation(firstCol, lastCol, k, firstRowW, firstColW, shift); - +#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE + ArrayXr tmp2 = (m_computed.block(firstCol+shift, firstCol+shift, n, n)).jacobiSvd().singularValues(); + std::cout << "\n\nj1 = " << tmp1.transpose().format(bdcsvdfmt) << "\n"; + std::cout << "j2 = " << tmp2.transpose().format(bdcsvdfmt) << "\n\n"; + std::cout << "err: " << ((tmp1-tmp2).abs()>1e-12*tmp2.abs()).transpose() << "\n"; + static int count = 0; + std::cout << "# " << ++count << "\n\n"; + assert((tmp1-tmp2).matrix().norm() < 1e-14*tmp2.matrix().norm()); +// assert(count<681); +// assert(((tmp1-tmp2).abs()<1e-13*tmp2.abs()).all()); +#endif + // Third part: compute SVD of combined matrix MatrixXr UofSVD, VofSVD; VectorType singVals; @@ -522,13 +541,24 @@ 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<actual_n;++k) + if(col0(k)!=0) + perm(m++) = k; + perm.conservativeResize(m); ArrayXr shifts(n), mus(n), zhat(n); @@ -539,12 +569,23 @@ #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"; std::cout << " sing-val: " << singVals.transpose() << "\n"; std::cout << " mu: " << mus.transpose() << "\n"; std::cout << " shift: " << shifts.transpose() << "\n"; + + { + Index actual_n = n; + while(actual_n>1 && col0(actual_n-1)==0) --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"; + 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"; + } #endif #ifdef EIGEN_BDCSVD_SANITY_CHECKS @@ -554,24 +595,16 @@ #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"; - { - Index actual_n = n; - while(actual_n>1 && col0(actual_n-1)==0) --actual_n; - std::cout << "\n\n mus: " << mus.head(actual_n).transpose() << "\n\n"; - std::cout << " check1: " << ((singVals.array()-(shifts+mus)) / singVals.array()).head(actual_n).transpose() << "\n\n"; - std::cout << " check2: " << ((singVals.array()-diag) / singVals.array()).head(actual_n).transpose() << "\n\n"; - std::cout << " check3: " << ((diag.segment(1,actual_n-1)-singVals.head(actual_n-1).array()) / singVals.head(actual_n-1).array()).transpose() << "\n\n\n"; - } #endif - + #ifdef EIGEN_BDCSVD_SANITY_CHECKS 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"; @@ -586,23 +619,48 @@ assert(m_computed.allFinite()); #endif + // Because of deflation, the singular values might not be completely sorted. + // Fortunately, reordering them is a O(n) problem + for(Index i=0; i<actual_n-1; ++i) + { + if(singVals(i)>singVals(i+1)) + { + using std::swap; + swap(singVals(i),singVals(i+1)); + U.col(i).swap(U.col(i+1)); + if(m_compV) V.col(i).swap(V.col(i+1)); + } + } + // Reverse order so that singular values in increased order // Because of deflation, the zeros singular-values are already at the end - Index actual_n = n; - while(actual_n>1 && singVals(actual_n-1)==0) --actual_n; singVals.head(actual_n).reverseInPlace(); U.leftCols(actual_n) = U.leftCols(actual_n).rowwise().reverse().eval(); // FIXME this requires a temporary if (m_compV) V.leftCols(actual_n) = V.leftCols(actual_n).rowwise().reverse().eval(); // FIXME this requires a temporary + +#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE + JacobiSVD<MatrixXr> jsvd(m_computed.block(firstCol, firstCol, n, n) ); + std::cout << " * j: " << jsvd.singularValues().transpose() << "\n\n"; + std::cout << " * sing-val: " << singVals.transpose() << "\n"; +// std::cout << " * err: " << ((jsvd.singularValues()-singVals)>1e-13*singVals.norm()).transpose() << "\n"; +#endif } template <typename MatrixType> -typename BDCSVD<MatrixType>::RealScalar BDCSVD<MatrixType>::secularEq(RealScalar mu, const ArrayXr& col0, const ArrayXr& diag, const ArrayXr& diagShifted, RealScalar shift, Index n) +typename BDCSVD<MatrixType>::RealScalar BDCSVD<MatrixType>::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<m; ++i) + { + Index j = perm(i); + res += numext::abs2(col0(j)) / ((diagShifted(j) - mu) * (diag(j) + shift + mu)); + } + return res; } template <typename MatrixType> -void BDCSVD<MatrixType>::computeSingVals(const ArrayXr& col0, const ArrayXr& diag, +void BDCSVD<MatrixType>::computeSingVals(const ArrayXr& col0, const ArrayXr& diag, const ArrayXi &perm, VectorType& singVals, ArrayXr& shifts, ArrayXr& mus) { using std::abs; @@ -612,26 +670,16 @@ Index n = col0.size(); Index actual_n = n; while(actual_n>1 && col0(actual_n-1)==0) --actual_n; -// Index m = 0; -// Array<Index,1,Dynamic> perm(actual_n); -// { -// for(Index k=0;k<actual_n;++k) -// { -// if(col0(k)!=0) -// perm(m++) = k; -// } -// } -// perm.conservativeResize(m); - + for (Index k = 0; k < n; ++k) { if (col0(k) == 0 || actual_n==1) { // if col0(k) == 0, then entry is deflated, so singular value is on diagonal // if actual_n==1, then the deflated problem is already diagonalized - singVals(k) = diag(k); + singVals(k) = k==0 ? col0(0) : diag(k); mus(k) = 0; - shifts(k) = diag(k); + shifts(k) = k==0 ? col0(0) : diag(k); continue; } @@ -650,8 +698,22 @@ // first decide whether it's closer to the left end or the right end RealScalar mid = left + (right-left) / 2; - RealScalar fMid = 1 + (col0.square() / ((diag + mid) * (diag - mid))).head(actual_n).sum(); - + RealScalar fMid = secularEq(mid, col0, diag, perm, diag, 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"; +#endif RealScalar shift = (k == actual_n-1 || fMid > 0) ? left : right; // measure everything relative to shift @@ -671,8 +733,8 @@ 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); @@ -682,20 +744,26 @@ // rational interpolation: fit a function of the form a / mu + b through the two previous // iterates and use its zero to compute the next iterate bool useBisection = fPrev*fCur>0; - while (abs(muCur - muPrev) > 8 * NumTraits<RealScalar>::epsilon() * (max)(abs(muCur), abs(muPrev)) && abs(fCur - fPrev)>NumTraits<RealScalar>::epsilon() && !useBisection) + while (fCur!=0 && abs(muCur - muPrev) > 8 * NumTraits<RealScalar>::epsilon() * (max)(abs(muCur), abs(muPrev)) && abs(fCur - fPrev)>NumTraits<RealScalar>::epsilon() && !useBisection) { ++m_numIters; + // Find a and b such that the function f(mu) = a / mu + b matches the current and previous samples. RealScalar a = (fCur - fPrev) / (1/muCur - 1/muPrev); RealScalar b = fCur - a / muCur; - + // And find mu such that f(mu)==0: + RealScalar muZero = -a/b; + RealScalar fZero = secularEq(muZero, col0, diag, perm, diagShifted, shift); + muPrev = muCur; fPrev = fCur; - muCur = -a / b; - fCur = secularEq(muCur, col0, diag, diagShifted, shift, actual_n); + muCur = muZero; + fCur = fZero; + if (shift == left && (muCur < 0 || muCur > right - left)) useBisection = true; if (shift == right && (muCur < -(right - left) || muCur > 0)) useBisection = true; + if (abs(fCur)>abs(fPrev)) useBisection = true; } // fall back on bisection method if rational interpolation did not work @@ -710,7 +778,7 @@ leftShifted = RealScalar(1)/NumTraits<RealScalar>::highest(); // I don't understand why the case k==0 would be special there: // if (k == 0) rightShifted = right - left; else - rightShifted = (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 } else { @@ -718,11 +786,11 @@ rightShifted = -RealScalar(1)/NumTraits<RealScalar>::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) + if(!(fLeft * fRight<0)) std::cout << k << " : " << fLeft << " * " << fRight << " == " << fLeft * fRight << " ; " << left << " - " << right << " -> " << leftShifted << " " << rightShifted << " shift=" << shift << "\n"; #endif eigen_internal_assert(fLeft * fRight < 0); @@ -730,7 +798,7 @@ while (rightShifted - leftShifted > 2 * NumTraits<RealScalar>::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; @@ -762,28 +830,13 @@ // zhat is perturbation of col0 for which singular vectors can be computed stably (see Section 3.1) template <typename MatrixType> void BDCSVD<MatrixType>::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<Index,1,Dynamic> perm(actual_n); - { - for(Index k=0;k<actual_n;++k) - { - if(col0(k)!=0) - perm(m++) = k; - } - } - perm.conservativeResize(m); + Index m = perm.size(); + Index last = perm(m-1); // The offset permits to skip deflated entries while computing zhat for (Index k = 0; k < n; ++k) @@ -794,7 +847,7 @@ { // 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<m; ++l) { @@ -805,12 +858,13 @@ prod *= ((singVals(j)+dk) / ((diag(i)+dk))) * ((mus(j)+(shifts(j)-dk)) / ((diag(i)-dk))); #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)) << ") / (" << (diag(i)+dk) << " * " << (diag(i)-dk) << ")\n"; + 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 } } #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; @@ -821,26 +875,11 @@ // compute singular vectors template <typename MatrixType> void BDCSVD<MatrixType>::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<Index,1,Dynamic> perm(actual_n); - { - for(Index k=0;k<actual_n;++k) - { - if(zhat(k)!=0) - perm(m++) = k; - } - } - perm.conservativeResize(m); + Index m = perm.size(); for (Index k = 0; k < n; ++k) { @@ -863,8 +902,6 @@ if (m_compV) { V.col(k).setZero(); -// for(Index i=1;i<actual_n;++i) -// V(i,k) = diag(i) * zhat(i) / (((diag(i) - shifts(k)) - mus(k)) )/( (diag(i) + singVals[k])); for(Index l=1;l<m;++l) { Index i = perm(l); @@ -908,7 +945,7 @@ // page 13 -// i,j >= 1, i != j and |di - dj| < epsilon * norm2(M) +// i,j >= 1, i!=j and |di - dj| < epsilon * norm2(M) // We apply two rotations to have zj = 0; // TODO deflation44 is still broken and not properly tested template <typename MatrixType> @@ -940,18 +977,13 @@ 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<RealScalar> J(c,s); - if (m_compU) - { - m_naiveU.middleRows(firstColu, size).applyOnTheRight(firstColu + i, firstColu + j, J); - } - if (m_compV) - { - m_naiveU.middleRows(firstRowW, size-1).applyOnTheRight(firstColW + i, firstColW + j, J.transpose()); - } + JacobiRotation<RealScalar> J(c,-s); + if (m_compU) m_naiveU.middleRows(firstColu, size+1).applyOnTheRight(firstColu + i, firstColu + j, J); + else m_naiveU.applyOnTheRight(firstColu+i, firstColu+j, J); + if (m_compV) m_naiveV.middleRows(firstRowW, size).applyOnTheRight(firstColW + i, firstColW + j, J); }// end deflation 44 @@ -964,43 +996,49 @@ using std::max; const Index length = lastCol + 1 - firstCol; + Block<MatrixXr,Dynamic,1> col0(m_computed, firstCol+shift, firstCol+shift, length, 1); + Diagonal<MatrixXr> fulldiag(m_computed); + VectorBlock<Diagonal<MatrixXr>,Dynamic> diag(fulldiag, firstCol+shift, length); + + RealScalar maxDiag = diag.tail((std::max)(Index(1),length-1)).cwiseAbs().maxCoeff(); + RealScalar epsilon_strict = NumTraits<RealScalar>::epsilon() * maxDiag; + RealScalar epsilon_coarse = 8 * NumTraits<RealScalar>::epsilon() * (max)(col0.cwiseAbs().maxCoeff(), maxDiag); + #ifdef EIGEN_BDCSVD_SANITY_CHECKS assert(m_naiveU.allFinite()); assert(m_naiveV.allFinite()); assert(m_computed.allFinite()); #endif - - Block<MatrixXr,Dynamic,1> col0(m_computed, firstCol+shift, firstCol+shift, length, 1); - Diagonal<MatrixXr> fulldiag(m_computed); - VectorBlock<Diagonal<MatrixXr>,Dynamic> diag(fulldiag, firstCol+shift, length); - - RealScalar epsilon = 8 * NumTraits<RealScalar>::epsilon() * (max)(col0.cwiseAbs().maxCoeff(), diag.cwiseAbs().maxCoeff()); + +#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE + std::cout << "\ndeflate:" << diag.head(k+1).transpose() << " | " << diag.segment(k+1,length-k-1).transpose() << "\n"; +#endif //condition 4.1 - if (diag(0) < epsilon) + if (diag(0) < epsilon_coarse) { #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE - std::cout << "deflation 4.1, because " << diag(0) << " < " << epsilon << "\n"; + std::cout << "deflation 4.1, because " << diag(0) << " < " << epsilon_coarse << "\n"; #endif - diag(0) = epsilon; + diag(0) = epsilon_coarse; } //condition 4.2 for (Index i=1;i<length;++i) - if (abs(col0(i)) < epsilon) + if (abs(col0(i)) < epsilon_strict) { #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE - std::cout << "deflation 4.2, set z(" << i << ") to zero because " << abs(col0(i)) << " < " << epsilon << " (diag(" << i << ")=" << diag(i) << ")\n"; + std::cout << "deflation 4.2, set z(" << i << ") to zero because " << abs(col0(i)) << " < " << epsilon_strict << " (diag(" << i << ")=" << diag(i) << ")\n"; #endif col0(i) = 0; } //condition 4.3 for (Index i=1;i<length; i++) - if (diag(i) < epsilon) + if (diag(i) < epsilon_coarse) { #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE - std::cout << "deflation 4.3, cancel z(" << i << ") because " << diag(i) << " < " << epsilon << " (z[" << i << "]=" << col0(i) << ")\n"; + std::cout << "deflation 4.3, cancel z(" << i << ")=" << col0(i) << " because diag(" << i << ")=" << diag(i) << " < " << epsilon_coarse << "\n"; #endif deflation43(firstCol, shift, i, length); } @@ -1010,14 +1048,22 @@ assert(m_naiveV.allFinite()); assert(m_computed.allFinite()); #endif - +#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE + std::cout << "to be sorted: " << diag.transpose() << "\n\n"; +#endif { + // 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 + bool total_deflation = (col0.tail(length-1).array()==RealScalar(0)).all(); + // 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. Index *permutation = new Index[length]; // FIXME avoid repeated dynamic memory allocation { permutation[0] = 0; Index p = 1; + + // Move deflated diagonal entries at the end. for(Index i=1; i<length; ++i) if(diag(i)==0) permutation[p++] = i; @@ -1032,6 +1078,22 @@ } } + // If we have a total deflation, then we have to insert diag(0) at the right place + if(total_deflation) + { + for(Index i=1; i<length; ++i) + { + Index pi = permutation[i]; + if(diag(pi)==0 || diag(0)<diag(pi)) + permutation[i-1] = permutation[i]; + else + { + permutation[i-1] = 0; + break; + } + } + } + // Current index of each col, and current column of each index Index *realInd = new Index[length]; // FIXME avoid repeated dynamic memory allocation Index *realCol = new Index[length]; // FIXME avoid repeated dynamic memory allocation @@ -1042,15 +1104,15 @@ realInd[pos] = pos; } - for(Index i = 1; i < length; i++) + for(Index i = total_deflation?0:1; i < length; i++) { - const Index pi = permutation[length - i]; + const Index pi = permutation[length - (total_deflation ? i+1 : i)]; const Index J = realCol[pi]; using std::swap; - // swap diaognal and first column entries: + // swap diagonal and first column entries: swap(diag(i), diag(J)); - swap(col0(i), col0(J)); + if(i!=0 && J!=0) swap(col0(i), col0(J)); // change columns if (m_compU) m_naiveU.col(firstCol+i).segment(firstCol, length + 1).swap(m_naiveU.col(firstCol+J).segment(firstCol, length + 1)); @@ -1068,23 +1130,31 @@ delete[] realInd; delete[] realCol; } +#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE + std::cout << "sorted: " << diag.transpose().format(bdcsvdfmt) << "\n"; + std::cout << " : " << col0.transpose() << "\n\n"; +#endif + + //condition 4.4 + { + 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<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"; +#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); + } + } #ifdef EIGEN_BDCSVD_SANITY_CHECKS - for(int k=2;k<length;++k) - assert(diag(k-1)<=diag(k) || diag(k)==0); + for(Index j=2;j<length;++j) + assert(diag(j-1)<=diag(j) || diag(j)==0); #endif - //condition 4.4 - for (Index i = 1; i+1<length && diag(i+1)!=0 && col0(i+1)!=0;i++) - if ((diag(i+1) - diag(i)) < NumTraits<RealScalar>::epsilon()*diag(i+1)) -// if ((diag(i+1) - diag(i)) < epsilon) - { -#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE - std::cout << "deflation 4.4 with i = " << i << " because " << (diag(i+1) - diag(i)) << " < " << epsilon << "\n"; -#endif - deflation44(firstCol, firstCol + shift, firstRowW, firstColW, i, i + 1, length); - } - #ifdef EIGEN_BDCSVD_SANITY_CHECKS assert(m_naiveU.allFinite()); assert(m_naiveV.allFinite());