20 #ifndef EIGEN_BDCSVD_H 21 #define EIGEN_BDCSVD_H 25 #ifdef EIGEN_BDCSVD_SANITY_CHECKS 26 #undef eigen_internal_assert 27 #define eigen_internal_assert(X) assert(X); 32 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE 33 IOFormat bdcsvdfmt(8, 0,
", ",
"\n",
" [",
"]");
36 template<
typename _MatrixType>
class BDCSVD;
40 template<
typename _MatrixType>
41 struct traits<
BDCSVD<_MatrixType> >
44 typedef _MatrixType MatrixType;
72 template<
typename _MatrixType>
83 typedef _MatrixType MatrixType;
84 typedef typename MatrixType::Scalar Scalar;
88 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
89 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
90 DiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime, ColsAtCompileTime),
91 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
92 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
93 MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(MaxRowsAtCompileTime, MaxColsAtCompileTime),
94 MatrixOptions = MatrixType::Options
99 typedef typename Base::SingularValuesType SingularValuesType;
114 BDCSVD() : m_algoswap(16), m_isTranspose(false), m_compU(false), m_compV(false), m_numIters(0)
125 : m_algoswap(16), m_numIters(0)
127 allocate(rows, cols, computationOptions);
140 BDCSVD(
const MatrixType& matrix,
unsigned int computationOptions = 0)
141 : m_algoswap(16), m_numIters(0)
143 compute(matrix, computationOptions);
160 BDCSVD& compute(
const MatrixType& matrix,
unsigned int computationOptions);
170 return compute(matrix, this->m_computationOptions);
173 void setSwitchSize(
int s)
175 eigen_assert(s>3 &&
"BDCSVD the size of the algo switch has to be greater than 3");
180 void allocate(
Index rows,
Index cols,
unsigned int computationOptions);
182 void computeSVDofM(
Index firstCol,
Index n, MatrixXr& U, VectorType& singVals, MatrixXr& V);
183 void computeSingVals(
const ArrayRef& col0,
const ArrayRef& diag,
const IndicesRef& perm, VectorType& singVals, ArrayRef shifts, ArrayRef mus);
184 void perturbCol0(
const ArrayRef& col0,
const ArrayRef& diag,
const IndicesRef& perm,
const VectorType& singVals,
const ArrayRef& shifts,
const ArrayRef& mus, ArrayRef zhat);
185 void computeSingVecs(
const ArrayRef& zhat,
const ArrayRef& diag,
const IndicesRef& perm,
const VectorType& singVals,
const ArrayRef& shifts,
const ArrayRef& mus, MatrixXr& U, MatrixXr& V);
189 template<
typename HouseholderU,
typename HouseholderV,
typename NaiveU,
typename NaiveV>
190 void copyUV(
const HouseholderU &householderU,
const HouseholderV &householderV,
const NaiveU &naiveU,
const NaiveV &naivev);
192 static RealScalar secularEq(RealScalar x,
const ArrayRef& col0,
const ArrayRef& diag,
const IndicesRef &perm,
const ArrayRef& diagShifted, RealScalar shift);
195 MatrixXr m_naiveU, m_naiveV;
199 ArrayXi m_workspaceI;
201 bool m_isTranspose, m_compU, m_compV;
203 using Base::m_singularValues;
204 using Base::m_diagSize;
205 using Base::m_computeFullU;
206 using Base::m_computeFullV;
207 using Base::m_computeThinU;
208 using Base::m_computeThinV;
209 using Base::m_matrixU;
210 using Base::m_matrixV;
212 using Base::m_isInitialized;
213 using Base::m_nonzeroSingularValues;
221 template<
typename MatrixType>
224 m_isTranspose = (cols > rows);
226 if (Base::allocate(rows, cols, computationOptions))
229 m_computed = MatrixXr::Zero(m_diagSize + 1, m_diagSize );
230 m_compU = computeV();
231 m_compV = computeU();
233 std::swap(m_compU, m_compV);
235 if (m_compU) m_naiveU = MatrixXr::Zero(m_diagSize + 1, m_diagSize + 1 );
236 else m_naiveU = MatrixXr::Zero(2, m_diagSize + 1 );
238 if (m_compV) m_naiveV = MatrixXr::Zero(m_diagSize, m_diagSize);
240 m_workspace.resize((m_diagSize+1)*(m_diagSize+1)*3);
241 m_workspaceI.resize(3*m_diagSize);
244 template<
typename MatrixType>
247 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE 248 std::cout <<
"\n\n\n======================================================================================================================\n\n\n";
250 allocate(matrix.rows(), matrix.cols(), computationOptions);
253 const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)();
256 if(matrix.cols() < m_algoswap)
260 m_isInitialized =
true;
261 m_info = jsvd.
info();
263 if(computeU()) m_matrixU = jsvd.
matrixU();
264 if(computeV()) m_matrixV = jsvd.
matrixV();
272 RealScalar scale = matrix.cwiseAbs().template maxCoeff<PropagateNaN>();
273 if (!(numext::isfinite)(scale)) {
274 m_isInitialized =
true;
279 if(scale==Literal(0)) scale = Literal(1);
281 if (m_isTranspose) copy = matrix.
adjoint()/scale;
282 else copy = matrix/scale;
286 internal::UpperBidiagonalization<MatrixX> bid(copy);
292 m_computed.topRows(m_diagSize) = bid.bidiagonal().toDenseMatrix().transpose();
293 m_computed.template bottomRows<1>().setZero();
294 divide(0, m_diagSize - 1, 0, 0, 0);
296 m_isInitialized =
true;
301 for (
int i=0; i<m_diagSize; i++)
303 RealScalar a =
abs(m_computed.coeff(i, i));
304 m_singularValues.coeffRef(i) = a * scale;
307 m_nonzeroSingularValues = i;
308 m_singularValues.tail(m_diagSize - i - 1).setZero();
311 else if (i == m_diagSize - 1)
313 m_nonzeroSingularValues = i + 1;
318 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE 322 if(m_isTranspose) copyUV(bid.householderV(), bid.householderU(), m_naiveV, m_naiveU);
323 else copyUV(bid.householderU(), bid.householderV(), m_naiveU, m_naiveV);
325 m_isInitialized =
true;
330 template<
typename MatrixType>
331 template<
typename HouseholderU,
typename HouseholderV,
typename NaiveU,
typename NaiveV>
332 void BDCSVD<MatrixType>::copyUV(
const HouseholderU &householderU,
const HouseholderV &householderV,
const NaiveU &naiveU,
const NaiveV &naiveV)
337 Index Ucols = m_computeThinU ? m_diagSize : householderU.cols();
339 m_matrixU.topLeftCorner(m_diagSize, m_diagSize) = naiveV.template cast<Scalar>().topLeftCorner(m_diagSize, m_diagSize);
340 householderU.applyThisOnTheLeft(m_matrixU);
344 Index Vcols = m_computeThinV ? m_diagSize : householderV.cols();
346 m_matrixV.topLeftCorner(m_diagSize, m_diagSize) = naiveU.template cast<Scalar>().topLeftCorner(m_diagSize, m_diagSize);
347 householderV.applyThisOnTheLeft(m_matrixV);
359 template<
typename MatrixType>
373 for(
Index j=0; j<n; ++j)
375 if( (A.col(j).head(n1).array()!=Literal(0)).any() )
377 A1.col(k1) = A.col(j).head(n1);
378 B1.row(k1) = B.row(j);
381 if( (A.col(j).tail(n2).array()!=Literal(0)).any() )
383 A2.col(k2) = A.col(j).tail(n2);
384 B2.row(k2) = B.row(j);
389 A.topRows(n1).noalias() = A1.leftCols(k1) * B1.topRows(k1);
390 A.bottomRows(n2).noalias() = A2.leftCols(k2) * B2.topRows(k2);
410 template<
typename MatrixType>
417 const Index n = lastCol - firstCol + 1;
419 const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)();
423 RealScalar lambda, phi, c0, s0;
434 m_naiveU.block(firstCol, firstCol, n + 1, n + 1).real() = b.
matrixU();
437 m_naiveU.row(0).segment(firstCol, n + 1).real() = b.
matrixU().row(0);
438 m_naiveU.row(1).segment(firstCol, n + 1).real() = b.
matrixU().row(n);
440 if (m_compV) m_naiveV.block(firstRowW, firstColW, n, n).real() = b.
matrixV();
441 m_computed.block(firstCol + shift, firstCol + shift, n + 1, n).setZero();
442 m_computed.diagonal().segment(firstCol + shift, n) = b.
singularValues().head(n);
446 alphaK = m_computed(firstCol + k, firstCol + k);
447 betaK = m_computed(firstCol + k + 1, firstCol + k);
451 divide(k + 1 + firstCol, lastCol, k + 1 + firstRowW, k + 1 + firstColW, shift);
453 divide(firstCol, k - 1 + firstCol, firstRowW, firstColW + 1, shift + 1);
458 lambda = m_naiveU(firstCol + k, firstCol + k);
459 phi = m_naiveU(firstCol + k + 1, lastCol + 1);
463 lambda = m_naiveU(1, firstCol + k);
464 phi = m_naiveU(0, lastCol + 1);
466 r0 =
sqrt((
abs(alphaK * lambda) *
abs(alphaK * lambda)) +
abs(betaK * phi) *
abs(betaK * phi));
469 l = m_naiveU.row(firstCol + k).segment(firstCol, k);
470 f = m_naiveU.row(firstCol + k + 1).segment(firstCol + k + 1, n - k - 1);
474 l = m_naiveU.row(1).segment(firstCol, k);
475 f = m_naiveU.row(0).segment(firstCol + k + 1, n - k - 1);
477 if (m_compV) m_naiveV(firstRowW+k, firstColW) = Literal(1);
485 c0 = alphaK * lambda / r0;
486 s0 = betaK * phi / r0;
489 #ifdef EIGEN_BDCSVD_SANITY_CHECKS 490 assert(m_naiveU.allFinite());
491 assert(m_naiveV.allFinite());
492 assert(m_computed.allFinite());
497 MatrixXr q1 (m_naiveU.col(firstCol + k).segment(firstCol, k + 1));
499 for (
Index i = firstCol + k - 1; i >= firstCol; i--)
500 m_naiveU.col(i + 1).segment(firstCol, k + 1) = m_naiveU.col(i).segment(firstCol, k + 1);
502 m_naiveU.col(firstCol).segment( firstCol, k + 1) = (q1 * c0);
504 m_naiveU.col(lastCol + 1).segment(firstCol, k + 1) = (q1 * ( - s0));
506 m_naiveU.col(firstCol).segment(firstCol + k + 1, n - k) = m_naiveU.col(lastCol + 1).segment(firstCol + k + 1, n - k) * s0;
508 m_naiveU.col(lastCol + 1).segment(firstCol + k + 1, n - k) *= c0;
512 RealScalar q1 = m_naiveU(0, firstCol + k);
514 for (
Index i = firstCol + k - 1; i >= firstCol; i--)
515 m_naiveU(0, i + 1) = m_naiveU(0, i);
517 m_naiveU(0, firstCol) = (q1 * c0);
519 m_naiveU(0, lastCol + 1) = (q1 * ( - s0));
521 m_naiveU(1, firstCol) = m_naiveU(1, lastCol + 1) *s0;
523 m_naiveU(1, lastCol + 1) *= c0;
524 m_naiveU.row(1).segment(firstCol + 1, k).setZero();
525 m_naiveU.row(0).segment(firstCol + k + 1, n - k - 1).setZero();
528 #ifdef EIGEN_BDCSVD_SANITY_CHECKS 529 assert(m_naiveU.allFinite());
530 assert(m_naiveV.allFinite());
531 assert(m_computed.allFinite());
534 m_computed(firstCol + shift, firstCol + shift) = r0;
535 m_computed.col(firstCol + shift).segment(firstCol + shift + 1, k) = alphaK * l.transpose().real();
536 m_computed.col(firstCol + shift).segment(firstCol + shift + k + 1, n - k - 1) = betaK * f.transpose().real();
538 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE 539 ArrayXr tmp1 = (m_computed.block(firstCol+shift, firstCol+shift, n, n)).jacobiSvd().singularValues();
542 deflation(firstCol, lastCol, k, firstRowW, firstColW, shift);
543 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE 544 ArrayXr tmp2 = (m_computed.block(firstCol+shift, firstCol+shift, n, n)).jacobiSvd().singularValues();
545 std::cout <<
"\n\nj1 = " << tmp1.transpose().format(bdcsvdfmt) <<
"\n";
546 std::cout <<
"j2 = " << tmp2.transpose().format(bdcsvdfmt) <<
"\n\n";
547 std::cout <<
"err: " << ((tmp1-tmp2).
abs()>1e-12*tmp2.abs()).transpose() <<
"\n";
548 static int count = 0;
549 std::cout <<
"# " << ++count <<
"\n\n";
550 assert((tmp1-tmp2).matrix().norm() < 1e-14*tmp2.matrix().norm());
558 computeSVDofM(firstCol + shift, n, UofSVD, singVals, VofSVD);
560 #ifdef EIGEN_BDCSVD_SANITY_CHECKS 561 assert(UofSVD.allFinite());
562 assert(VofSVD.allFinite());
566 structured_update(m_naiveU.block(firstCol, firstCol, n + 1, n + 1), UofSVD, (n+2)/2);
570 tmp.noalias() = m_naiveU.middleCols(firstCol, n+1) * UofSVD;
571 m_naiveU.middleCols(firstCol, n + 1) = tmp;
574 if (m_compV) structured_update(m_naiveV.block(firstRowW, firstColW, n, n), VofSVD, (n+1)/2);
576 #ifdef EIGEN_BDCSVD_SANITY_CHECKS 577 assert(m_naiveU.allFinite());
578 assert(m_naiveV.allFinite());
579 assert(m_computed.allFinite());
582 m_computed.block(firstCol + shift, firstCol + shift, n, n).setZero();
583 m_computed.block(firstCol + shift, firstCol + shift, n, n).diagonal() = singVals;
594 template <
typename MatrixType>
597 const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)();
599 ArrayRef col0 = m_computed.col(firstCol).segment(firstCol, n);
600 m_workspace.head(n) = m_computed.block(firstCol, firstCol, n, n).diagonal();
601 ArrayRef diag = m_workspace.head(n);
602 diag(0) = Literal(0);
607 if (m_compV) V.
resize(n, n);
609 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE 610 if (col0.hasNaN() || diag.hasNaN())
611 std::cout <<
"\n\nHAS NAN\n\n";
618 while(actual_n>1 && diag(actual_n-1)==Literal(0)) {--actual_n; eigen_internal_assert(col0(actual_n)==Literal(0)); }
620 for(
Index k=0;k<actual_n;++k)
621 if(
abs(col0(k))>considerZero)
622 m_workspaceI(m++) = k;
629 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE 630 std::cout <<
"computeSVDofM using:\n";
631 std::cout <<
" z: " << col0.transpose() <<
"\n";
632 std::cout <<
" d: " << diag.transpose() <<
"\n";
636 computeSingVals(col0, diag, perm, singVals, shifts, mus);
638 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE 639 std::cout <<
" j: " << (m_computed.block(firstCol, firstCol, n, n)).jacobiSvd().singularValues().transpose().reverse() <<
"\n\n";
640 std::cout <<
" sing-val: " << singVals.transpose() <<
"\n";
641 std::cout <<
" mu: " << mus.transpose() <<
"\n";
642 std::cout <<
" shift: " << shifts.transpose() <<
"\n";
645 std::cout <<
"\n\n mus: " << mus.head(actual_n).transpose() <<
"\n\n";
646 std::cout <<
" check1 (expect0) : " << ((singVals.array()-(shifts+mus)) / singVals.array()).head(actual_n).transpose() <<
"\n\n";
647 assert((((singVals.array()-(shifts+mus)) / singVals.array()).head(actual_n) >= 0).
all());
648 std::cout <<
" check2 (>0) : " << ((singVals.array()-diag) / singVals.array()).head(actual_n).transpose() <<
"\n\n";
649 assert((((singVals.array()-diag) / singVals.array()).head(actual_n) >= 0).
all());
653 #ifdef EIGEN_BDCSVD_SANITY_CHECKS 654 assert(singVals.allFinite());
655 assert(mus.allFinite());
656 assert(shifts.allFinite());
660 perturbCol0(col0, diag, perm, singVals, shifts, mus, zhat);
661 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE 662 std::cout <<
" zhat: " << zhat.transpose() <<
"\n";
665 #ifdef EIGEN_BDCSVD_SANITY_CHECKS 666 assert(zhat.allFinite());
669 computeSingVecs(zhat, diag, perm, singVals, shifts, mus, U, V);
671 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE 672 std::cout <<
"U^T U: " << (U.transpose() * U -
MatrixXr(MatrixXr::Identity(U.cols(),U.cols()))).norm() <<
"\n";
673 std::cout <<
"V^T V: " << (V.transpose() * V -
MatrixXr(MatrixXr::Identity(V.cols(),V.cols()))).norm() <<
"\n";
676 #ifdef EIGEN_BDCSVD_SANITY_CHECKS 677 assert(m_naiveU.allFinite());
678 assert(m_naiveV.allFinite());
679 assert(m_computed.allFinite());
680 assert(U.allFinite());
681 assert(V.allFinite());
688 for(
Index i=0; i<actual_n-1; ++i)
690 if(singVals(i)>singVals(i+1))
693 swap(singVals(i),singVals(i+1));
694 U.col(i).swap(U.col(i+1));
695 if(m_compV) V.col(i).swap(V.col(i+1));
699 #ifdef EIGEN_BDCSVD_SANITY_CHECKS 701 bool singular_values_sorted = (((singVals.segment(1,actual_n-1)-singVals.head(actual_n-1))).array() >= 0).
all();
702 if(!singular_values_sorted)
703 std::cout <<
"Singular values are not sorted: " << singVals.segment(1,actual_n).transpose() <<
"\n";
704 assert(singular_values_sorted);
710 singVals.head(actual_n).reverseInPlace();
711 U.leftCols(actual_n).rowwise().reverseInPlace();
712 if (m_compV) V.leftCols(actual_n).rowwise().reverseInPlace();
714 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE 716 std::cout <<
" * j: " << jsvd.singularValues().transpose() <<
"\n\n";
717 std::cout <<
" * sing-val: " << singVals.transpose() <<
"\n";
722 template <
typename MatrixType>
725 Index m = perm.size();
726 RealScalar res = Literal(1);
727 for(
Index i=0; i<m; ++i)
732 res += (col0(j) / (diagShifted(j) - mu)) * (col0(j) / (diag(j) + shift + mu));
738 template <
typename MatrixType>
746 Index n = col0.size();
750 while(actual_n>1 && col0(actual_n-1)==Literal(0)) --actual_n;
752 for (
Index k = 0; k < n; ++k)
754 if (col0(k) == Literal(0) || actual_n==1)
758 singVals(k) = k==0 ? col0(0) : diag(k);
760 shifts(k) = k==0 ? col0(0) : diag(k);
765 RealScalar left = diag(k);
768 right = (diag(actual_n-1) + col0.matrix().norm());
775 while(col0(l)==Literal(0)) { ++l; eigen_internal_assert(l<actual_n); }
780 RealScalar mid = left + (right-left) / Literal(2);
781 RealScalar fMid = secularEq(mid, col0, diag, perm, diag, Literal(0));
782 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE 783 std::cout <<
"right-left = " << right-left <<
"\n";
786 std::cout <<
" = " << secularEq(left+RealScalar(0.000001)*(right-left), col0, diag, perm, diag, 0)
787 <<
" " << secularEq(left+RealScalar(0.1) *(right-left), col0, diag, perm, diag, 0)
788 <<
" " << secularEq(left+RealScalar(0.2) *(right-left), col0, diag, perm, diag, 0)
789 <<
" " << secularEq(left+RealScalar(0.3) *(right-left), col0, diag, perm, diag, 0)
790 <<
" " << secularEq(left+RealScalar(0.4) *(right-left), col0, diag, perm, diag, 0)
791 <<
" " << secularEq(left+RealScalar(0.49) *(right-left), col0, diag, perm, diag, 0)
792 <<
" " << secularEq(left+RealScalar(0.5) *(right-left), col0, diag, perm, diag, 0)
793 <<
" " << secularEq(left+RealScalar(0.51) *(right-left), col0, diag, perm, diag, 0)
794 <<
" " << secularEq(left+RealScalar(0.6) *(right-left), col0, diag, perm, diag, 0)
795 <<
" " << secularEq(left+RealScalar(0.7) *(right-left), col0, diag, perm, diag, 0)
796 <<
" " << secularEq(left+RealScalar(0.8) *(right-left), col0, diag, perm, diag, 0)
797 <<
" " << secularEq(left+RealScalar(0.9) *(right-left), col0, diag, perm, diag, 0)
798 <<
" " << secularEq(left+RealScalar(0.999999)*(right-left), col0, diag, perm, diag, 0) <<
"\n";
800 RealScalar shift = (k == actual_n-1 || fMid > Literal(0)) ? left : right;
804 diagShifted = diag - shift;
809 RealScalar midShifted = (right - left) / RealScalar(2);
811 midShifted = -midShifted;
812 RealScalar fMidShifted = secularEq(midShifted, col0, diag, perm, diagShifted, shift);
816 shift = fMidShifted > Literal(0) ? left : right;
817 diagShifted = diag - shift;
822 RealScalar muPrev, muCur;
825 muPrev = (right - left) * RealScalar(0.1);
826 if (k == actual_n-1) muCur = right - left;
827 else muCur = (right - left) * RealScalar(0.5);
831 muPrev = -(right - left) * RealScalar(0.1);
832 muCur = -(right - left) * RealScalar(0.5);
835 RealScalar fPrev = secularEq(muPrev, col0, diag, perm, diagShifted, shift);
836 RealScalar fCur = secularEq(muCur, col0, diag, perm, diagShifted, shift);
837 if (
abs(fPrev) <
abs(fCur))
845 bool useBisection = fPrev*fCur>Literal(0);
851 RealScalar a = (fCur - fPrev) / (Literal(1)/muCur - Literal(1)/muPrev);
852 RealScalar b = fCur - a / muCur;
854 RealScalar muZero = -a/b;
855 RealScalar fZero = secularEq(muZero, col0, diag, perm, diagShifted, shift);
857 #ifdef EIGEN_BDCSVD_SANITY_CHECKS 858 assert((numext::isfinite)(fZero));
866 if (shift == left && (muCur < Literal(0) || muCur > right - left)) useBisection =
true;
867 if (shift == right && (muCur < -(right - left) || muCur > Literal(0))) useBisection =
true;
868 if (
abs(fCur)>
abs(fPrev)) useBisection =
true;
874 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE 875 std::cout <<
"useBisection for k = " << k <<
", actual_n = " << actual_n <<
"\n";
877 RealScalar leftShifted, rightShifted;
882 leftShifted = numext::maxi<RealScalar>( (std::numeric_limits<RealScalar>::min)(), Literal(2) *
abs(col0(k)) /
sqrt((std::numeric_limits<RealScalar>::max)()) );
885 eigen_internal_assert( (numext::isfinite)( (col0(k)/leftShifted)*(col0(k)/(diag(k)+shift+leftShifted)) ) );
888 rightShifted = (k==actual_n-1) ? right : ((right - left) * RealScalar(0.51));
892 leftShifted = -(right - left) * RealScalar(0.51);
894 rightShifted = -numext::maxi<RealScalar>( (std::numeric_limits<RealScalar>::min)(),
abs(col0(k+1)) /
sqrt((std::numeric_limits<RealScalar>::max)()) );
896 rightShifted = -(std::numeric_limits<RealScalar>::min)();
899 RealScalar fLeft = secularEq(leftShifted, col0, diag, perm, diagShifted, shift);
900 eigen_internal_assert(fLeft<Literal(0));
902 #if defined EIGEN_INTERNAL_DEBUGGING || defined EIGEN_BDCSVD_SANITY_CHECKS 903 RealScalar fRight = secularEq(rightShifted, col0, diag, perm, diagShifted, shift);
906 #ifdef EIGEN_BDCSVD_SANITY_CHECKS 907 if(!(numext::isfinite)(fLeft))
908 std::cout <<
"f(" << leftShifted <<
") =" << fLeft <<
" ; " << left <<
" " << shift <<
" " << right <<
"\n";
909 assert((numext::isfinite)(fLeft));
911 if(!(numext::isfinite)(fRight))
912 std::cout <<
"f(" << rightShifted <<
") =" << fRight <<
" ; " << left <<
" " << shift <<
" " << right <<
"\n";
916 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE 917 if(!(fLeft * fRight<0))
919 std::cout <<
"f(leftShifted) using leftShifted=" << leftShifted <<
" ; diagShifted(1:10):" << diagShifted.head(10).transpose() <<
"\n ; " 920 <<
"left==shift=" << bool(left==shift) <<
" ; left-shift = " << (left-shift) <<
"\n";
921 std::cout <<
"k=" << k <<
", " << fLeft <<
" * " << fRight <<
" == " << fLeft * fRight <<
" ; " 922 <<
"[" << left <<
" .. " << right <<
"] -> [" << leftShifted <<
" " << rightShifted <<
"], shift=" << shift
923 <<
" , f(right)=" << secularEq(0, col0, diag, perm, diagShifted, shift)
924 <<
" == " << secularEq(right, col0, diag, perm, diag, 0) <<
" == " << fRight <<
"\n";
927 eigen_internal_assert(fLeft * fRight < Literal(0));
933 RealScalar midShifted = (leftShifted + rightShifted) / Literal(2);
934 fMid = secularEq(midShifted, col0, diag, perm, diagShifted, shift);
935 eigen_internal_assert((numext::isfinite)(fMid));
937 if (fLeft * fMid < Literal(0))
939 rightShifted = midShifted;
943 leftShifted = midShifted;
947 muCur = (leftShifted + rightShifted) / Literal(2);
955 muCur = (right - left) * RealScalar(0.5);
961 singVals[k] = shift + muCur;
965 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE 967 std::cout <<
"found " << singVals[k] <<
" == " << shift <<
" + " << muCur <<
" from " << diag(k) <<
" .. " << diag(k+1) <<
"\n";
969 #ifdef EIGEN_BDCSVD_SANITY_CHECKS 970 assert(k==0 || singVals[k]>=singVals[k-1]);
971 assert(singVals[k]>=diag(k));
984 template <
typename MatrixType>
990 Index n = col0.size();
991 Index m = perm.size();
997 Index lastIdx = perm(m-1);
999 for (
Index k = 0; k < n; ++k)
1001 if (col0(k) == Literal(0))
1002 zhat(k) = Literal(0);
1006 RealScalar dk = diag(k);
1007 RealScalar prod = (singVals(lastIdx) + dk) * (mus(lastIdx) + (shifts(lastIdx) - dk));
1008 #ifdef EIGEN_BDCSVD_SANITY_CHECKS 1010 std::cout <<
"k = " << k <<
" ; z(k)=" << col0(k) <<
", diag(k)=" << dk <<
"\n";
1011 std::cout <<
"prod = " <<
"(" << singVals(lastIdx) <<
" + " << dk <<
") * (" << mus(lastIdx) <<
" + (" << shifts(lastIdx) <<
" - " << dk <<
"))" <<
"\n";
1012 std::cout <<
" = " << singVals(lastIdx) + dk <<
" * " << mus(lastIdx) + (shifts(lastIdx) - dk) <<
"\n";
1017 for(
Index l = 0; l<m; ++l)
1022 #ifdef EIGEN_BDCSVD_SANITY_CHECKS 1023 if(i>=k && (l==0 || l-1>=m))
1025 std::cout <<
"Error in perturbCol0\n";
1026 std::cout <<
" " << k <<
"/" << n <<
" " << l <<
"/" << m <<
" " << i <<
"/" << n <<
" ; " << col0(k) <<
" " << diag(k) <<
" " <<
"\n";
1027 std::cout <<
" " <<diag(i) <<
"\n";
1028 Index j = (i<k ) ? i : perm(l-1);
1029 std::cout <<
" " <<
"j=" << j <<
"\n";
1032 Index j = i<k ? i : perm(l-1);
1033 #ifdef EIGEN_BDCSVD_SANITY_CHECKS 1034 if(!(dk!=Literal(0) || diag(i)!=Literal(0)))
1036 std::cout <<
"k=" << k <<
", i=" << i <<
", l=" << l <<
", perm.size()=" << perm.size() <<
"\n";
1038 assert(dk!=Literal(0) || diag(i)!=Literal(0));
1040 prod *= ((singVals(j)+dk) / ((diag(i)+dk))) * ((mus(j)+(shifts(j)-dk)) / ((diag(i)-dk)));
1041 #ifdef EIGEN_BDCSVD_SANITY_CHECKS 1044 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE 1045 if(i!=k && numext::abs(((singVals(j)+dk)*(mus(j)+(shifts(j)-dk)))/((diag(i)+dk)*(diag(i)-dk)) - 1) > 0.9 )
1046 std::cout <<
" " << ((singVals(j)+dk)*(mus(j)+(shifts(j)-dk)))/((diag(i)+dk)*(diag(i)-dk)) <<
" == (" << (singVals(j)+dk) <<
" * " << (mus(j)+(shifts(j)-dk))
1047 <<
") / (" << (diag(i)+dk) <<
" * " << (diag(i)-dk) <<
")\n";
1051 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE 1052 std::cout <<
"zhat(" << k <<
") = sqrt( " << prod <<
") ; " << (singVals(lastIdx) + dk) <<
" * " << mus(lastIdx) + shifts(lastIdx) <<
" - " << dk <<
"\n";
1054 RealScalar tmp =
sqrt(prod);
1055 #ifdef EIGEN_BDCSVD_SANITY_CHECKS 1056 assert((numext::isfinite)(tmp));
1058 zhat(k) = col0(k) > Literal(0) ? RealScalar(tmp) : RealScalar(-tmp);
1064 template <
typename MatrixType>
1069 Index n = zhat.size();
1070 Index m = perm.size();
1072 for (
Index k = 0; k < n; ++k)
1074 if (zhat(k) == Literal(0))
1076 U.col(k) = VectorType::Unit(n+1, k);
1077 if (m_compV) V.col(k) = VectorType::Unit(n, k);
1082 for(
Index l=0;l<m;++l)
1085 U(i,k) = zhat(i)/(((diag(i) - shifts(k)) - mus(k)) )/( (diag(i) + singVals[k]));
1087 U(n,k) = Literal(0);
1088 U.col(k).normalize();
1093 for(
Index l=1;l<m;++l)
1096 V(i,k) = diag(i) * zhat(i) / (((diag(i) - shifts(k)) - mus(k)) )/( (diag(i) + singVals[k]));
1098 V(0,k) = Literal(-1);
1099 V.col(k).normalize();
1103 U.col(n) = VectorType::Unit(n+1, n);
1110 template <
typename MatrixType>
1116 Index start = firstCol + shift;
1117 RealScalar c = m_computed(start, start);
1118 RealScalar s = m_computed(start+i, start);
1119 RealScalar r = numext::hypot(c,s);
1120 if (r == Literal(0))
1122 m_computed(start+i, start+i) = Literal(0);
1125 m_computed(start,start) = r;
1126 m_computed(start+i, start) = Literal(0);
1127 m_computed(start+i, start+i) = Literal(0);
1130 if (m_compU) m_naiveU.middleRows(firstCol, size+1).applyOnTheRight(firstCol, firstCol+i, J);
1131 else m_naiveU.applyOnTheRight(firstCol, firstCol+i, J);
1139 template <
typename MatrixType>
1146 RealScalar c = m_computed(firstColm+i, firstColm);
1147 RealScalar s = m_computed(firstColm+j, firstColm);
1148 RealScalar r =
sqrt(numext::abs2(c) + numext::abs2(s));
1149 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE 1150 std::cout <<
"deflation 4.4: " << i <<
"," << j <<
" -> " << c <<
" " << s <<
" " << r <<
" ; " 1151 << m_computed(firstColm + i-1, firstColm) <<
" " 1152 << m_computed(firstColm + i, firstColm) <<
" " 1153 << m_computed(firstColm + i+1, firstColm) <<
" " 1154 << m_computed(firstColm + i+2, firstColm) <<
"\n";
1155 std::cout << m_computed(firstColm + i-1, firstColm + i-1) <<
" " 1156 << m_computed(firstColm + i, firstColm+i) <<
" " 1157 << m_computed(firstColm + i+1, firstColm+i+1) <<
" " 1158 << m_computed(firstColm + i+2, firstColm+i+2) <<
"\n";
1162 m_computed(firstColm + i, firstColm + i) = m_computed(firstColm + j, firstColm + j);
1167 m_computed(firstColm + i, firstColm) = r;
1168 m_computed(firstColm + j, firstColm + j) = m_computed(firstColm + i, firstColm + i);
1169 m_computed(firstColm + j, firstColm) = Literal(0);
1172 if (m_compU) m_naiveU.middleRows(firstColu, size+1).applyOnTheRight(firstColu + i, firstColu + j, J);
1173 else m_naiveU.applyOnTheRight(firstColu+i, firstColu+j, J);
1174 if (m_compV) m_naiveV.middleRows(firstRowW, size).applyOnTheRight(firstColW + i, firstColW + j, J);
1179 template <
typename MatrixType>
1184 const Index length = lastCol + 1 - firstCol;
1190 const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)();
1191 RealScalar maxDiag = diag.tail((std::max)(
Index(1),length-1)).cwiseAbs().maxCoeff();
1195 #ifdef EIGEN_BDCSVD_SANITY_CHECKS 1196 assert(m_naiveU.allFinite());
1197 assert(m_naiveV.allFinite());
1198 assert(m_computed.allFinite());
1201 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE 1202 std::cout <<
"\ndeflate:" << diag.head(k+1).transpose() <<
" | " << diag.segment(k+1,length-k-1).transpose() <<
"\n";
1206 if (diag(0) < epsilon_coarse)
1208 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE 1209 std::cout <<
"deflation 4.1, because " << diag(0) <<
" < " << epsilon_coarse <<
"\n";
1211 diag(0) = epsilon_coarse;
1215 for (
Index i=1;i<length;++i)
1216 if (
abs(col0(i)) < epsilon_strict)
1218 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE 1219 std::cout <<
"deflation 4.2, set z(" << i <<
") to zero because " <<
abs(col0(i)) <<
" < " << epsilon_strict <<
" (diag(" << i <<
")=" << diag(i) <<
")\n";
1221 col0(i) = Literal(0);
1225 for (
Index i=1;i<length; i++)
1226 if (diag(i) < epsilon_coarse)
1228 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE 1229 std::cout <<
"deflation 4.3, cancel z(" << i <<
")=" << col0(i) <<
" because diag(" << i <<
")=" << diag(i) <<
" < " << epsilon_coarse <<
"\n";
1231 deflation43(firstCol, shift, i, length);
1234 #ifdef EIGEN_BDCSVD_SANITY_CHECKS 1235 assert(m_naiveU.allFinite());
1236 assert(m_naiveV.allFinite());
1237 assert(m_computed.allFinite());
1239 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE 1240 std::cout <<
"to be sorted: " << diag.transpose() <<
"\n\n";
1241 std::cout <<
" : " << col0.transpose() <<
"\n\n";
1246 bool total_deflation = (col0.tail(length-1).array()<considerZero).
all();
1250 Index *permutation = m_workspaceI.data();
1256 for(
Index i=1; i<length; ++i)
1257 if(
abs(diag(i))<considerZero)
1258 permutation[p++] = i;
1261 for( ; p < length; ++p)
1263 if (i > k) permutation[p] = j++;
1264 else if (j >= length) permutation[p] = i++;
1265 else if (diag(i) < diag(j)) permutation[p] = j++;
1266 else permutation[p] = i++;
1273 for(
Index i=1; i<length; ++i)
1275 Index pi = permutation[i];
1276 if(
abs(diag(pi))<considerZero || diag(0)<diag(pi))
1277 permutation[i-1] = permutation[i];
1280 permutation[i-1] = 0;
1287 Index *realInd = m_workspaceI.data()+length;
1288 Index *realCol = m_workspaceI.data()+2*length;
1290 for(
int pos = 0; pos< length; pos++)
1296 for(
Index i = total_deflation?0:1; i < length; i++)
1298 const Index pi = permutation[length - (total_deflation ? i+1 : i)];
1299 const Index J = realCol[pi];
1303 swap(diag(i), diag(J));
1304 if(i!=0 && J!=0) swap(col0(i), col0(J));
1307 if (m_compU) m_naiveU.col(firstCol+i).segment(firstCol, length + 1).swap(m_naiveU.col(firstCol+J).segment(firstCol, length + 1));
1308 else m_naiveU.col(firstCol+i).segment(0, 2) .swap(m_naiveU.col(firstCol+J).segment(0, 2));
1309 if (m_compV) m_naiveV.col(firstColW + i).segment(firstRowW, length).swap(m_naiveV.col(firstColW + J).segment(firstRowW, length));
1312 const Index realI = realInd[i];
1319 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE 1320 std::cout <<
"sorted: " << diag.transpose().format(bdcsvdfmt) <<
"\n";
1321 std::cout <<
" : " << col0.transpose() <<
"\n\n";
1327 while(i>0 && (
abs(diag(i))<considerZero ||
abs(col0(i))<considerZero)) --i;
1331 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE 1332 std::cout <<
"deflation 4.4 with i = " << i <<
" because " << diag(i) <<
" - " << diag(i-1) <<
" == " << (diag(i) - diag(i-1)) <<
" < " <<
NumTraits<RealScalar>::epsilon()*maxDiag <<
"\n";
1334 eigen_internal_assert(
abs(diag(i) - diag(i-1))<epsilon_coarse &&
" diagonal entries are not properly sorted");
1335 deflation44(firstCol, firstCol + shift, firstRowW, firstColW, i-1, i, length);
1339 #ifdef EIGEN_BDCSVD_SANITY_CHECKS 1340 for(
Index j=2;j<length;++j)
1341 assert(diag(j-1)<=diag(j) ||
abs(diag(j))<considerZero);
1344 #ifdef EIGEN_BDCSVD_SANITY_CHECKS 1345 assert(m_naiveU.allFinite());
1346 assert(m_naiveV.allFinite());
1347 assert(m_computed.allFinite());
1357 template<
typename Derived>
BDCSVD & compute(const MatrixType &matrix)
Method performing the decomposition of given matrix using current options.
Definition: BDCSVD.h:168
Definition: Constants.h:393
BDCSVD< PlainObject > bdcSvd(unsigned int computationOptions=0) const
Definition: BDCSVD.h:1359
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SVDBase.h:236
static const Eigen::internal::all_t all
Definition: IndexedViewHelper.h:171
const AdjointReturnType adjoint() const
Definition: Transpose.h:221
Matrix< Type, Dynamic, Dynamic > MatrixX
[c++11]
Definition: Matrix.h:543
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(const Eigen::ArrayBase< Derived > &x)
const MatrixUType & matrixU() const
Definition: SVDBase.h:101
Namespace containing all symbols from the Eigen library.
Definition: Core:141
Rotation given by a cosine-sine pair.
Definition: Jacobi.h:34
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:232
void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:271
Derived & setZero(Index size)
Definition: CwiseNullaryOp.h:562
Expression of a fixed-size or dynamic-size sub-vector.
Definition: VectorBlock.h:56
BDCSVD(Index rows, Index cols, unsigned int computationOptions=0)
Default Constructor with memory preallocation.
Definition: BDCSVD.h:124
Base class of SVD algorithms.
Definition: SVDBase.h:62
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Definition: Constants.h:240
BDCSVD & compute(const MatrixType &matrix, unsigned int computationOptions)
Method performing the decomposition of given matrix using custom options.
Definition: BDCSVD.h:245
BDCSVD(const MatrixType &matrix, unsigned int computationOptions=0)
Constructor performing the decomposition of given matrix.
Definition: BDCSVD.h:140
Definition: Constants.h:449
class Bidiagonal Divide and Conquer SVD
Definition: BDCSVD.h:36
Definition: Constants.h:442
A matrix or vector expression mapping an existing expression.
Definition: Ref.h:281
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
Definition: Eigen_Colamd.h:50
static const IdentityReturnType Identity()
Definition: CwiseNullaryOp.h:799
Expression of a fixed-size or dynamic-size block.
Definition: Block.h:103
Index nonzeroSingularValues() const
Definition: SVDBase.h:136
BDCSVD()
Default Constructor.
Definition: BDCSVD.h:114
Two-sided Jacobi SVD decomposition of a rectangular matrix.
Definition: JacobiSVD.h:488
const SingularValuesType & singularValues() const
Definition: SVDBase.h:129
const MatrixVType & matrixV() const
Definition: SVDBase.h:117
Expression of a diagonal/subdiagonal/superdiagonal in a matrix.
Definition: Diagonal.h:63
const int Dynamic
Definition: Constants.h:22
Definition: Constants.h:397
Definition: Constants.h:446
Array< int, Dynamic, 1 > ArrayXi
Definition: Array.h:352