GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-11 11:19:02 +08:00
Code Issues Packages Projects Releases Wiki Activity

Fix total deflation issue in BDCSVD, when & only when M is already diagonal.

Browse Source
This commit is contained in:
Xinle Liu 2021-11-02 16:53:55 +00:00 committed by Rasmus Munk Larsen
parent 8f8c2ba2fe
commit 478a1bdda6
2 changed files with 158 additions and 121 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

228
Eigen/src/SVD/BDCSVD.h
View File

@ -1,9 +1,9 @@
// This file is part of Eigen, a lightweight C++ template library // This file is part of Eigen, a lightweight C++ template library
// for linear algebra. // for linear algebra.
// //
// We used the "A Divide-And-Conquer Algorithm for the Bidiagonal SVD" // We used the "A Divide-And-Conquer Algorithm for the Bidiagonal SVD"
// research report written by Ming Gu and Stanley C.Eisenstat // research report written by Ming Gu and Stanley C.Eisenstat
// The code variable names correspond to the names they used in their // The code variable names correspond to the names they used in their
// report // report
// //
// Copyright (C) 2013 Gauthier Brun <brun.gauthier@gmail.com> // Copyright (C) 2013 Gauthier Brun <brun.gauthier@gmail.com>
@ -29,12 +29,16 @@
#include "./InternalHeaderCheck.h" #include "./InternalHeaderCheck.h"
#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
#include <iostream>
#endif
namespace Eigen { namespace Eigen {
#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
IOFormat bdcsvdfmt(8, 0, ", ", "\n", " [", "]"); IOFormat bdcsvdfmt(8, 0, ", ", "\n", " [", "]");
#endif #endif
template<typename MatrixType_> class BDCSVD; template<typename MatrixType_> class BDCSVD;
namespace internal { namespace internal {
@ -44,11 +48,11 @@ struct traits<BDCSVD<MatrixType_> >
: traits<MatrixType_> : traits<MatrixType_>
{ {
typedef MatrixType_ MatrixType; typedef MatrixType_ MatrixType;
}; };
} // end namespace internal } // end namespace internal
/** \ingroup SVD_Module /** \ingroup SVD_Module
* *
* *
@ -75,31 +79,31 @@ template<typename MatrixType_>
class BDCSVD : public SVDBase<BDCSVD<MatrixType_> > class BDCSVD : public SVDBase<BDCSVD<MatrixType_> >
{ {
typedef SVDBase<BDCSVD> Base; typedef SVDBase<BDCSVD> Base;
public: public:
using Base::rows; using Base::rows;
using Base::cols; using Base::cols;
using Base::computeU; using Base::computeU;
using Base::computeV; using Base::computeV;
typedef MatrixType_ MatrixType; typedef MatrixType_ MatrixType;
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
typedef typename NumTraits<RealScalar>::Literal Literal; typedef typename NumTraits<RealScalar>::Literal Literal;
enum { enum {
RowsAtCompileTime = MatrixType::RowsAtCompileTime, RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ColsAtCompileTime = MatrixType::ColsAtCompileTime, ColsAtCompileTime = MatrixType::ColsAtCompileTime,
DiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime, ColsAtCompileTime), DiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime, ColsAtCompileTime),
MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(MaxRowsAtCompileTime, MaxColsAtCompileTime), MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(MaxRowsAtCompileTime, MaxColsAtCompileTime),
MatrixOptions = MatrixType::Options MatrixOptions = MatrixType::Options
}; };
typedef typename Base::MatrixUType MatrixUType; typedef typename Base::MatrixUType MatrixUType;
typedef typename Base::MatrixVType MatrixVType; typedef typename Base::MatrixVType MatrixVType;
typedef typename Base::SingularValuesType SingularValuesType; typedef typename Base::SingularValuesType SingularValuesType;
typedef Matrix<Scalar, Dynamic, Dynamic, ColMajor> MatrixX; typedef Matrix<Scalar, Dynamic, Dynamic, ColMajor> MatrixX;
typedef Matrix<RealScalar, Dynamic, Dynamic, ColMajor> MatrixXr; typedef Matrix<RealScalar, Dynamic, Dynamic, ColMajor> MatrixXr;
typedef Matrix<RealScalar, Dynamic, 1> VectorType; typedef Matrix<RealScalar, Dynamic, 1> VectorType;
@ -133,7 +137,7 @@ public:
* *
* \param matrix the matrix to decompose * \param matrix the matrix to decompose
* \param computationOptions optional parameter allowing to specify if you want full or thin U or V unitaries to be computed. * \param computationOptions optional parameter allowing to specify if you want full or thin U or V unitaries to be computed.
* By default, none is computed. This is a bit - field, the possible bits are #ComputeFullU, #ComputeThinU, * By default, none is computed. This is a bit - field, the possible bits are #ComputeFullU, #ComputeThinU,
* #ComputeFullV, #ComputeThinV. * #ComputeFullV, #ComputeThinV.
* *
* Thin unitaries are only available if your matrix type has a Dynamic number of columns (for example MatrixXf). They also are not * Thin unitaries are only available if your matrix type has a Dynamic number of columns (for example MatrixXf). They also are not
@ -145,15 +149,15 @@ public:
compute(matrix, computationOptions); compute(matrix, computationOptions);
} }
~BDCSVD() ~BDCSVD()
{ {
} }
/** \brief Method performing the decomposition of given matrix using custom options. /** \brief Method performing the decomposition of given matrix using custom options.
* *
* \param matrix the matrix to decompose * \param matrix the matrix to decompose
* \param computationOptions optional parameter allowing to specify if you want full or thin U or V unitaries to be computed. * \param computationOptions optional parameter allowing to specify if you want full or thin U or V unitaries to be computed.
* By default, none is computed. This is a bit - field, the possible bits are #ComputeFullU, #ComputeThinU, * By default, none is computed. This is a bit - field, the possible bits are #ComputeFullU, #ComputeThinU,
* #ComputeFullV, #ComputeThinV. * #ComputeFullV, #ComputeThinV.
* *
* Thin unitaries are only available if your matrix type has a Dynamic number of columns (for example MatrixXf). They also are not * Thin unitaries are only available if your matrix type has a Dynamic number of columns (for example MatrixXf). They also are not
@ -172,12 +176,12 @@ public:
return compute(matrix, this->m_computationOptions); return compute(matrix, this->m_computationOptions);
} }
void setSwitchSize(int s) void setSwitchSize(int s)
{ {
eigen_assert(s>3 && "BDCSVD the size of the algo switch has to be greater than 3"); eigen_assert(s>=3 && "BDCSVD the size of the algo switch has to be at least 3.");
m_algoswap = s; m_algoswap = s;
} }
private: private:
void allocate(Index rows, Index cols, unsigned int computationOptions); void allocate(Index rows, Index cols, unsigned int computationOptions);
void divide(Index firstCol, Index lastCol, Index firstRowW, Index firstColW, Index shift); void divide(Index firstCol, Index lastCol, Index firstRowW, Index firstColW, Index shift);
@ -201,7 +205,7 @@ protected:
ArrayXi m_workspaceI; ArrayXi m_workspaceI;
int m_algoswap; int m_algoswap;
bool m_isTranspose, m_compU, m_compV; bool m_isTranspose, m_compU, m_compV;
using Base::m_singularValues; using Base::m_singularValues;
using Base::m_diagSize; using Base::m_diagSize;
using Base::m_computeFullU; using Base::m_computeFullU;
@ -214,7 +218,7 @@ protected:
using Base::m_isInitialized; using Base::m_isInitialized;
using Base::m_nonzeroSingularValues; using Base::m_nonzeroSingularValues;
public: public:
int m_numIters; int m_numIters;
}; //end class BDCSVD }; //end class BDCSVD
@ -227,24 +231,24 @@ void BDCSVD<MatrixType>::allocate(Eigen::Index rows, Eigen::Index cols, unsigned
if (Base::allocate(rows, cols, computationOptions)) if (Base::allocate(rows, cols, computationOptions))
return; return;
m_computed = MatrixXr::Zero(m_diagSize + 1, m_diagSize ); m_computed = MatrixXr::Zero(m_diagSize + 1, m_diagSize );
m_compU = computeV(); m_compU = computeV();
m_compV = computeU(); m_compV = computeU();
if (m_isTranspose) if (m_isTranspose)
std::swap(m_compU, m_compV); std::swap(m_compU, m_compV);
if (m_compU) m_naiveU = MatrixXr::Zero(m_diagSize + 1, m_diagSize + 1 ); if (m_compU) m_naiveU = MatrixXr::Zero(m_diagSize + 1, m_diagSize + 1 );
else m_naiveU = MatrixXr::Zero(2, m_diagSize + 1 ); else m_naiveU = MatrixXr::Zero(2, m_diagSize + 1 );
if (m_compV) m_naiveV = MatrixXr::Zero(m_diagSize, m_diagSize); if (m_compV) m_naiveV = MatrixXr::Zero(m_diagSize, m_diagSize);
m_workspace.resize((m_diagSize+1)*(m_diagSize+1)*3); m_workspace.resize((m_diagSize+1)*(m_diagSize+1)*3);
m_workspaceI.resize(3*m_diagSize); m_workspaceI.resize(3*m_diagSize);
}// end allocate }// end allocate
template<typename MatrixType> template<typename MatrixType>
BDCSVD<MatrixType>& BDCSVD<MatrixType>::compute(const MatrixType& matrix, unsigned int computationOptions) BDCSVD<MatrixType>& BDCSVD<MatrixType>::compute(const MatrixType& matrix, unsigned int computationOptions)
{ {
#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
std::cout << "\n\n\n======================================================================================================================\n\n\n"; std::cout << "\n\n\n======================================================================================================================\n\n\n";
@ -253,7 +257,7 @@ BDCSVD<MatrixType>& BDCSVD<MatrixType>::compute(const MatrixType& matrix, unsign
using std::abs; using std::abs;
const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)(); const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)();
//**** step -1 - If the problem is too small, directly falls back to JacobiSVD and return //**** step -1 - If the problem is too small, directly falls back to JacobiSVD and return
if(matrix.cols() < m_algoswap) if(matrix.cols() < m_algoswap)
{ {
@ -269,7 +273,7 @@ BDCSVD<MatrixType>& BDCSVD<MatrixType>::compute(const MatrixType& matrix, unsign
} }
return *this; return *this;
} }
//**** step 0 - Copy the input matrix and apply scaling to reduce over/under-flows //**** step 0 - Copy the input matrix and apply scaling to reduce over/under-flows
RealScalar scale = matrix.cwiseAbs().template maxCoeff<PropagateNaN>(); RealScalar scale = matrix.cwiseAbs().template maxCoeff<PropagateNaN>();
if (!(numext::isfinite)(scale)) { if (!(numext::isfinite)(scale)) {
@ -282,7 +286,7 @@ BDCSVD<MatrixType>& BDCSVD<MatrixType>::compute(const MatrixType& matrix, unsign
MatrixX copy; MatrixX copy;
if (m_isTranspose) copy = matrix.adjoint()/scale; if (m_isTranspose) copy = matrix.adjoint()/scale;
else copy = matrix/scale; else copy = matrix/scale;
//**** step 1 - Bidiagonalization //**** step 1 - Bidiagonalization
// FIXME this line involves temporaries // FIXME this line involves temporaries
internal::UpperBidiagonalization<MatrixX> bid(copy); internal::UpperBidiagonalization<MatrixX> bid(copy);
@ -298,7 +302,7 @@ BDCSVD<MatrixType>& BDCSVD<MatrixType>::compute(const MatrixType& matrix, unsign
m_isInitialized = true; m_isInitialized = true;
return *this; return *this;
} }
//**** step 3 - Copy singular values and vectors //**** step 3 - Copy singular values and vectors
for (int i=0; i<m_diagSize; i++) for (int i=0; i<m_diagSize; i++)
{ {
@ -387,7 +391,7 @@ void BDCSVD<MatrixType>::structured_update(Block<MatrixXr,Dynamic,Dynamic> A, co
++k2; ++k2;
} }
} }
A.topRows(n1).noalias() = A1.leftCols(k1) * B1.topRows(k1); A.topRows(n1).noalias() = A1.leftCols(k1) * B1.topRows(k1);
A.bottomRows(n2).noalias() = A2.leftCols(k2) * B2.topRows(k2); A.bottomRows(n2).noalias() = A2.leftCols(k2) * B2.topRows(k2);
} }
@ -399,15 +403,15 @@ void BDCSVD<MatrixType>::structured_update(Block<MatrixXr,Dynamic,Dynamic> A, co
} }
} }
// The divide algorithm is done "in place", we are always working on subsets of the same matrix. The divide methods takes as argument the // The divide algorithm is done "in place", we are always working on subsets of the same matrix. The divide methods takes as argument the
// place of the submatrix we are currently working on. // place of the submatrix we are currently working on.
//@param firstCol : The Index of the first column of the submatrix of m_computed and for m_naiveU; //@param firstCol : The Index of the first column of the submatrix of m_computed and for m_naiveU;
//@param lastCol : The Index of the last column of the submatrix of m_computed and for m_naiveU; //@param lastCol : The Index of the last column of the submatrix of m_computed and for m_naiveU;
// lastCol + 1 - firstCol is the size of the submatrix. // lastCol + 1 - firstCol is the size of the submatrix.
//@param firstRowW : The Index of the first row of the matrix W that we are to change. (see the reference paper section 1 for more information on W) //@param firstRowW : The Index of the first row of the matrix W that we are to change. (see the reference paper section 1 for more information on W)
//@param firstRowW : Same as firstRowW with the column. //@param firstColW : Same as firstRowW with the column.
//@param shift : Each time one takes the left submatrix, one must add 1 to the shift. Why? Because! We actually want the last column of the U submatrix //@param shift : Each time one takes the left submatrix, one must add 1 to the shift. Why? Because! We actually want the last column of the U submatrix
// to become the first column (*coeff) and to shift all the other columns to the right. There are more details on the reference paper. // to become the first column (*coeff) and to shift all the other columns to the right. There are more details on the reference paper.
template<typename MatrixType> template<typename MatrixType>
void BDCSVD<MatrixType>::divide(Eigen::Index firstCol, Eigen::Index lastCol, Eigen::Index firstRowW, Eigen::Index firstColW, Eigen::Index shift) void BDCSVD<MatrixType>::divide(Eigen::Index firstCol, Eigen::Index lastCol, Eigen::Index firstRowW, Eigen::Index firstColW, Eigen::Index shift)
@ -420,11 +424,11 @@ void BDCSVD<MatrixType>::divide(Eigen::Index firstCol, Eigen::Index lastCol, Eig
const Index k = n/2; const Index k = n/2;
const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)(); const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)();
RealScalar alphaK; RealScalar alphaK;
RealScalar betaK; RealScalar betaK;
RealScalar r0; RealScalar r0;
RealScalar lambda, phi, c0, s0; RealScalar lambda, phi, c0, s0;
VectorType l, f; VectorType l, f;
// We use the other algorithm which is more efficient for small // We use the other algorithm which is more efficient for small
// matrices. // matrices.
if (n < m_algoswap) if (n < m_algoswap)
{ {
@ -434,7 +438,7 @@ void BDCSVD<MatrixType>::divide(Eigen::Index firstCol, Eigen::Index lastCol, Eig
if (m_info != Success && m_info != NoConvergence) return; if (m_info != Success && m_info != NoConvergence) return;
if (m_compU) if (m_compU)
m_naiveU.block(firstCol, firstCol, n + 1, n + 1).real() = b.matrixU(); m_naiveU.block(firstCol, firstCol, n + 1, n + 1).real() = b.matrixU();
else else
{ {
m_naiveU.row(0).segment(firstCol, n + 1).real() = b.matrixU().row(0); m_naiveU.row(0).segment(firstCol, n + 1).real() = b.matrixU().row(0);
m_naiveU.row(1).segment(firstCol, n + 1).real() = b.matrixU().row(n); m_naiveU.row(1).segment(firstCol, n + 1).real() = b.matrixU().row(n);
@ -448,8 +452,8 @@ void BDCSVD<MatrixType>::divide(Eigen::Index firstCol, Eigen::Index lastCol, Eig
alphaK = m_computed(firstCol + k, firstCol + k); alphaK = m_computed(firstCol + k, firstCol + k);
betaK = m_computed(firstCol + k + 1, firstCol + k); betaK = m_computed(firstCol + k + 1, firstCol + k);
// The divide must be done in that order in order to have good results. Divide change the data inside the submatrices // The divide must be done in that order in order to have good results. Divide change the data inside the submatrices
// and the divide of the right submatrice reads one column of the left submatrice. That's why we need to treat the // and the divide of the right submatrice reads one column of the left submatrice. That's why we need to treat the
// right submatrix before the left one. // right submatrix before the left one.
divide(k + 1 + firstCol, lastCol, k + 1 + firstRowW, k + 1 + firstColW, shift); divide(k + 1 + firstCol, lastCol, k + 1 + firstRowW, k + 1 + firstColW, shift);
if (m_info != Success && m_info != NoConvergence) return; if (m_info != Success && m_info != NoConvergence) return;
divide(firstCol, k - 1 + firstCol, firstRowW, firstColW + 1, shift + 1); divide(firstCol, k - 1 + firstCol, firstRowW, firstColW + 1, shift + 1);
@ -459,8 +463,8 @@ void BDCSVD<MatrixType>::divide(Eigen::Index firstCol, Eigen::Index lastCol, Eig
{ {
lambda = m_naiveU(firstCol + k, firstCol + k); lambda = m_naiveU(firstCol + k, firstCol + k);
phi = m_naiveU(firstCol + k + 1, lastCol + 1); phi = m_naiveU(firstCol + k + 1, lastCol + 1);
} }
else else
{ {
lambda = m_naiveU(1, firstCol + k); lambda = m_naiveU(1, firstCol + k);
phi = m_naiveU(0, lastCol + 1); phi = m_naiveU(0, lastCol + 1);
@ -470,8 +474,8 @@ void BDCSVD<MatrixType>::divide(Eigen::Index firstCol, Eigen::Index lastCol, Eig
{ {
l = m_naiveU.row(firstCol + k).segment(firstCol, k); l = m_naiveU.row(firstCol + k).segment(firstCol, k);
f = m_naiveU.row(firstCol + k + 1).segment(firstCol + k + 1, n - k - 1); f = m_naiveU.row(firstCol + k + 1).segment(firstCol + k + 1, n - k - 1);
} }
else else
{ {
l = m_naiveU.row(1).segment(firstCol, k); l = m_naiveU.row(1).segment(firstCol, k);
f = m_naiveU.row(0).segment(firstCol + k + 1, n - k - 1); f = m_naiveU.row(0).segment(firstCol + k + 1, n - k - 1);
@ -487,52 +491,52 @@ void BDCSVD<MatrixType>::divide(Eigen::Index firstCol, Eigen::Index lastCol, Eig
c0 = alphaK * lambda / r0; c0 = alphaK * lambda / r0;
s0 = betaK * phi / r0; s0 = betaK * phi / r0;
} }
#ifdef EIGEN_BDCSVD_SANITY_CHECKS #ifdef EIGEN_BDCSVD_SANITY_CHECKS
assert(m_naiveU.allFinite()); assert(m_naiveU.allFinite());
assert(m_naiveV.allFinite()); assert(m_naiveV.allFinite());
assert(m_computed.allFinite()); assert(m_computed.allFinite());
#endif #endif
if (m_compU) if (m_compU)
{ {
MatrixXr q1 (m_naiveU.col(firstCol + k).segment(firstCol, k + 1)); MatrixXr q1 (m_naiveU.col(firstCol + k).segment(firstCol, k + 1));
// we shiftW Q1 to the right // we shiftW Q1 to the right
for (Index i = firstCol + k - 1; i >= firstCol; i--) for (Index i = firstCol + k - 1; i >= firstCol; i--)
m_naiveU.col(i + 1).segment(firstCol, k + 1) = m_naiveU.col(i).segment(firstCol, k + 1); m_naiveU.col(i + 1).segment(firstCol, k + 1) = m_naiveU.col(i).segment(firstCol, k + 1);
// we shift q1 at the left with a factor c0 // we shift q1 at the left with a factor c0
m_naiveU.col(firstCol).segment( firstCol, k + 1) = (q1 * c0); m_naiveU.col(firstCol).segment( firstCol, k + 1) = (q1 * c0);
// last column = q1 * - s0 // last column = q1 * - s0
m_naiveU.col(lastCol + 1).segment(firstCol, k + 1) = (q1 * ( - s0)); m_naiveU.col(lastCol + 1).segment(firstCol, k + 1) = (q1 * ( - s0));
// first column = q2 * s0 // first column = q2 * s0
m_naiveU.col(firstCol).segment(firstCol + k + 1, n - k) = m_naiveU.col(lastCol + 1).segment(firstCol + k + 1, n - k) * s0; m_naiveU.col(firstCol).segment(firstCol + k + 1, n - k) = m_naiveU.col(lastCol + 1).segment(firstCol + k + 1, n - k) * s0;
// q2 *= c0 // q2 *= c0
m_naiveU.col(lastCol + 1).segment(firstCol + k + 1, n - k) *= c0; m_naiveU.col(lastCol + 1).segment(firstCol + k + 1, n - k) *= c0;
} }
else else
{ {
RealScalar q1 = m_naiveU(0, firstCol + k); RealScalar q1 = m_naiveU(0, firstCol + k);
// we shift Q1 to the right // we shift Q1 to the right
for (Index i = firstCol + k - 1; i >= firstCol; i--) for (Index i = firstCol + k - 1; i >= firstCol; i--)
m_naiveU(0, i + 1) = m_naiveU(0, i); m_naiveU(0, i + 1) = m_naiveU(0, i);
// we shift q1 at the left with a factor c0 // we shift q1 at the left with a factor c0
m_naiveU(0, firstCol) = (q1 * c0); m_naiveU(0, firstCol) = (q1 * c0);
// last column = q1 * - s0 // last column = q1 * - s0
m_naiveU(0, lastCol + 1) = (q1 * ( - s0)); m_naiveU(0, lastCol + 1) = (q1 * ( - s0));
// first column = q2 * s0 // first column = q2 * s0
m_naiveU(1, firstCol) = m_naiveU(1, lastCol + 1) *s0; m_naiveU(1, firstCol) = m_naiveU(1, lastCol + 1) *s0;
// q2 *= c0 // q2 *= c0
m_naiveU(1, lastCol + 1) *= c0; m_naiveU(1, lastCol + 1) *= c0;
m_naiveU.row(1).segment(firstCol + 1, k).setZero(); m_naiveU.row(1).segment(firstCol + 1, k).setZero();
m_naiveU.row(0).segment(firstCol + k + 1, n - k - 1).setZero(); m_naiveU.row(0).segment(firstCol + k + 1, n - k - 1).setZero();
} }
#ifdef EIGEN_BDCSVD_SANITY_CHECKS #ifdef EIGEN_BDCSVD_SANITY_CHECKS
assert(m_naiveU.allFinite()); assert(m_naiveU.allFinite());
assert(m_naiveV.allFinite()); assert(m_naiveV.allFinite());
assert(m_computed.allFinite()); assert(m_computed.allFinite());
#endif #endif
m_computed(firstCol + shift, firstCol + shift) = r0; m_computed(firstCol + shift, firstCol + shift) = r0;
m_computed.col(firstCol + shift).segment(firstCol + shift + 1, k) = alphaK * l.transpose().real(); 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(); m_computed.col(firstCol + shift).segment(firstCol + shift + k + 1, n - k - 1) = betaK * f.transpose().real();
@ -553,17 +557,17 @@ void BDCSVD<MatrixType>::divide(Eigen::Index firstCol, Eigen::Index lastCol, Eig
// assert(count<681); // assert(count<681);
// assert(((tmp1-tmp2).abs()<1e-13*tmp2.abs()).all()); // assert(((tmp1-tmp2).abs()<1e-13*tmp2.abs()).all());
#endif #endif
// Third part: compute SVD of combined matrix // Third part: compute SVD of combined matrix
MatrixXr UofSVD, VofSVD; MatrixXr UofSVD, VofSVD;
VectorType singVals; VectorType singVals;
computeSVDofM(firstCol + shift, n, UofSVD, singVals, VofSVD); computeSVDofM(firstCol + shift, n, UofSVD, singVals, VofSVD);
#ifdef EIGEN_BDCSVD_SANITY_CHECKS #ifdef EIGEN_BDCSVD_SANITY_CHECKS
assert(UofSVD.allFinite()); assert(UofSVD.allFinite());
assert(VofSVD.allFinite()); assert(VofSVD.allFinite());
#endif #endif
if (m_compU) if (m_compU)
structured_update(m_naiveU.block(firstCol, firstCol, n + 1, n + 1), UofSVD, (n+2)/2); structured_update(m_naiveU.block(firstCol, firstCol, n + 1, n + 1), UofSVD, (n+2)/2);
else else
@ -572,15 +576,15 @@ void BDCSVD<MatrixType>::divide(Eigen::Index firstCol, Eigen::Index lastCol, Eig
tmp.noalias() = m_naiveU.middleCols(firstCol, n+1) * UofSVD; tmp.noalias() = m_naiveU.middleCols(firstCol, n+1) * UofSVD;
m_naiveU.middleCols(firstCol, n + 1) = tmp; m_naiveU.middleCols(firstCol, n + 1) = tmp;
} }
if (m_compV) structured_update(m_naiveV.block(firstRowW, firstColW, n, n), VofSVD, (n+1)/2); if (m_compV) structured_update(m_naiveV.block(firstRowW, firstColW, n, n), VofSVD, (n+1)/2);
#ifdef EIGEN_BDCSVD_SANITY_CHECKS #ifdef EIGEN_BDCSVD_SANITY_CHECKS
assert(m_naiveU.allFinite()); assert(m_naiveU.allFinite());
assert(m_naiveV.allFinite()); assert(m_naiveV.allFinite());
assert(m_computed.allFinite()); assert(m_computed.allFinite());
#endif #endif
m_computed.block(firstCol + shift, firstCol + shift, n, n).setZero(); m_computed.block(firstCol + shift, firstCol + shift, n, n).setZero();
m_computed.block(firstCol + shift, firstCol + shift, n, n).diagonal() = singVals; m_computed.block(firstCol + shift, firstCol + shift, n, n).diagonal() = singVals;
}// end divide }// end divide
@ -612,7 +616,7 @@ void BDCSVD<MatrixType>::computeSVDofM(Eigen::Index firstCol, Eigen::Index n, Ma
if (col0.hasNaN() || diag.hasNaN()) if (col0.hasNaN() || diag.hasNaN())
std::cout << "\n\nHAS NAN\n\n"; std::cout << "\n\nHAS NAN\n\n";
#endif #endif
// Many singular values might have been deflated, the zero ones have been moved to the end, // 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. // but others are interleaved and we must ignore them at this stage.
// To this end, let's compute a permutation skipping them: // To this end, let's compute a permutation skipping them:
@ -623,7 +627,7 @@ void BDCSVD<MatrixType>::computeSVDofM(Eigen::Index firstCol, Eigen::Index n, Ma
if(abs(col0(k))>considerZero) if(abs(col0(k))>considerZero)
m_workspaceI(m++) = k; m_workspaceI(m++) = k;
Map<ArrayXi> perm(m_workspaceI.data(),m); Map<ArrayXi> perm(m_workspaceI.data(),m);
Map<ArrayXr> shifts(m_workspace.data()+1*n, n); Map<ArrayXr> shifts(m_workspace.data()+1*n, n);
Map<ArrayXr> mus(m_workspace.data()+2*n, n); Map<ArrayXr> mus(m_workspace.data()+2*n, n);
Map<ArrayXr> zhat(m_workspace.data()+3*n, n); Map<ArrayXr> zhat(m_workspace.data()+3*n, n);
@ -633,16 +637,16 @@ void BDCSVD<MatrixType>::computeSVDofM(Eigen::Index firstCol, Eigen::Index n, Ma
std::cout << " z: " << col0.transpose() << "\n"; std::cout << " z: " << col0.transpose() << "\n";
std::cout << " d: " << diag.transpose() << "\n"; std::cout << " d: " << diag.transpose() << "\n";
#endif #endif
// Compute singVals, shifts, and mus // Compute singVals, shifts, and mus
computeSingVals(col0, diag, perm, singVals, shifts, mus); computeSingVals(col0, diag, perm, singVals, shifts, mus);
#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
std::cout << " j: " << (m_computed.block(firstCol, firstCol, n, n)).jacobiSvd().singularValues().transpose().reverse() << "\n\n"; 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 << " sing-val: " << singVals.transpose() << "\n";
std::cout << " mu: " << mus.transpose() << "\n"; std::cout << " mu: " << mus.transpose() << "\n";
std::cout << " shift: " << shifts.transpose() << "\n"; std::cout << " shift: " << shifts.transpose() << "\n";
{ {
std::cout << "\n\n mus: " << mus.head(actual_n).transpose() << "\n\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 << " check1 (expect0) : " << ((singVals.array()-(shifts+mus)) / singVals.array()).head(actual_n).transpose() << "\n\n";
@ -651,30 +655,30 @@ void BDCSVD<MatrixType>::computeSVDofM(Eigen::Index firstCol, Eigen::Index n, Ma
assert((((singVals.array()-diag) / singVals.array()).head(actual_n) >= 0).all()); assert((((singVals.array()-diag) / singVals.array()).head(actual_n) >= 0).all());
} }
#endif #endif
#ifdef EIGEN_BDCSVD_SANITY_CHECKS #ifdef EIGEN_BDCSVD_SANITY_CHECKS
assert(singVals.allFinite()); assert(singVals.allFinite());
assert(mus.allFinite()); assert(mus.allFinite());
assert(shifts.allFinite()); assert(shifts.allFinite());
#endif #endif
// Compute zhat // Compute zhat
perturbCol0(col0, diag, perm, singVals, shifts, mus, zhat); perturbCol0(col0, diag, perm, singVals, shifts, mus, zhat);
#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
std::cout << " zhat: " << zhat.transpose() << "\n"; std::cout << " zhat: " << zhat.transpose() << "\n";
#endif #endif
#ifdef EIGEN_BDCSVD_SANITY_CHECKS #ifdef EIGEN_BDCSVD_SANITY_CHECKS
assert(zhat.allFinite()); assert(zhat.allFinite());
#endif #endif
computeSingVecs(zhat, diag, perm, singVals, shifts, mus, U, V); computeSingVecs(zhat, diag, perm, singVals, shifts, mus, U, V);
#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
std::cout << "U^T U: " << (U.transpose() * U - MatrixXr(MatrixXr::Identity(U.cols(),U.cols()))).norm() << "\n"; std::cout << "U^T U: " << (U.transpose() * U - MatrixXr(MatrixXr::Identity(U.cols(),U.cols()))).norm() << "\n";
std::cout << "V^T V: " << (V.transpose() * V - MatrixXr(MatrixXr::Identity(V.cols(),V.cols()))).norm() << "\n"; std::cout << "V^T V: " << (V.transpose() * V - MatrixXr(MatrixXr::Identity(V.cols(),V.cols()))).norm() << "\n";
#endif #endif
#ifdef EIGEN_BDCSVD_SANITY_CHECKS #ifdef EIGEN_BDCSVD_SANITY_CHECKS
assert(m_naiveU.allFinite()); assert(m_naiveU.allFinite());
assert(m_naiveV.allFinite()); assert(m_naiveV.allFinite());
@ -684,7 +688,7 @@ void BDCSVD<MatrixType>::computeSVDofM(Eigen::Index firstCol, Eigen::Index n, Ma
// assert((U.transpose() * U - MatrixXr(MatrixXr::Identity(U.cols(),U.cols()))).norm() < 100*NumTraits<RealScalar>::epsilon() * n); // 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); // assert((V.transpose() * V - MatrixXr(MatrixXr::Identity(V.cols(),V.cols()))).norm() < 100*NumTraits<RealScalar>::epsilon() * n);
#endif #endif
// Because of deflation, the singular values might not be completely sorted. // Because of deflation, the singular values might not be completely sorted.
// Fortunately, reordering them is a O(n) problem // Fortunately, reordering them is a O(n) problem
for(Index i=0; i<actual_n-1; ++i) for(Index i=0; i<actual_n-1; ++i)
@ -706,13 +710,13 @@ void BDCSVD<MatrixType>::computeSVDofM(Eigen::Index firstCol, Eigen::Index n, Ma
assert(singular_values_sorted); assert(singular_values_sorted);
} }
#endif #endif
// Reverse order so that singular values in increased order // Reverse order so that singular values in increased order
// Because of deflation, the zeros singular-values are already at the end // Because of deflation, the zeros singular-values are already at the end
singVals.head(actual_n).reverseInPlace(); singVals.head(actual_n).reverseInPlace();
U.leftCols(actual_n).rowwise().reverseInPlace(); U.leftCols(actual_n).rowwise().reverseInPlace();
if (m_compV) V.leftCols(actual_n).rowwise().reverseInPlace(); if (m_compV) V.leftCols(actual_n).rowwise().reverseInPlace();
#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
JacobiSVD<MatrixXr> jsvd(m_computed.block(firstCol, firstCol, n, n) ); JacobiSVD<MatrixXr> jsvd(m_computed.block(firstCol, firstCol, n, n) );
std::cout << " * j: " << jsvd.singularValues().transpose() << "\n\n"; std::cout << " * j: " << jsvd.singularValues().transpose() << "\n\n";
@ -761,7 +765,7 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
mus(k) = Literal(0); mus(k) = Literal(0);
shifts(k) = k==0 ? col0(0) : diag(k); shifts(k) = k==0 ? col0(0) : diag(k);
continue; continue;
} }
// otherwise, use secular equation to find singular value // otherwise, use secular equation to find singular value
RealScalar left = diag(k); RealScalar left = diag(k);
@ -800,7 +804,7 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
<< " " << secularEq(left+RealScalar(0.999999)*(right-left), col0, diag, perm, diag, 0) << "\n"; << " " << secularEq(left+RealScalar(0.999999)*(right-left), col0, diag, perm, diag, 0) << "\n";
#endif #endif
RealScalar shift = (k == actual_n-1 || fMid > Literal(0)) ? left : right; RealScalar shift = (k == actual_n-1 || fMid > Literal(0)) ? left : right;
// measure everything relative to shift // measure everything relative to shift
Map<ArrayXr> diagShifted(m_workspace.data()+4*n, n); Map<ArrayXr> diagShifted(m_workspace.data()+4*n, n);
diagShifted = diag - shift; diagShifted = diag - shift;
@ -819,7 +823,7 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
diagShifted = diag - shift; diagShifted = diag - shift;
} }
} }
// initial guess // initial guess
RealScalar muPrev, muCur; RealScalar muPrev, muCur;
if (shift == left) if (shift == left)
@ -859,12 +863,12 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
#ifdef EIGEN_BDCSVD_SANITY_CHECKS #ifdef EIGEN_BDCSVD_SANITY_CHECKS
assert((numext::isfinite)(fZero)); assert((numext::isfinite)(fZero));
#endif #endif
muPrev = muCur; muPrev = muCur;
fPrev = fCur; fPrev = fCur;
muCur = muZero; muCur = muZero;
fCur = fZero; fCur = fZero;
if (shift == left && (muCur < Literal(0) || muCur > right - left)) useBisection = true; if (shift == left && (muCur < Literal(0) || muCur > right - left)) useBisection = true;
if (shift == right && (muCur < -(right - left) || muCur > Literal(0))) useBisection = true; if (shift == right && (muCur < -(right - left) || muCur > Literal(0))) useBisection = true;
if (abs(fCur)>abs(fPrev)) useBisection = true; if (abs(fCur)>abs(fPrev)) useBisection = true;
@ -901,7 +905,7 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
RealScalar fLeft = secularEq(leftShifted, col0, diag, perm, diagShifted, shift); RealScalar fLeft = secularEq(leftShifted, col0, diag, perm, diagShifted, shift);
eigen_internal_assert(fLeft<Literal(0)); eigen_internal_assert(fLeft<Literal(0));
#if defined EIGEN_INTERNAL_DEBUGGING || defined EIGEN_BDCSVD_SANITY_CHECKS #if defined EIGEN_BDCSVD_DEBUG_VERBOSE || defined EIGEN_BDCSVD_SANITY_CHECKS || defined EIGEN_INTERNAL_DEBUGGING
RealScalar fRight = secularEq(rightShifted, col0, diag, perm, diagShifted, shift); RealScalar fRight = secularEq(rightShifted, col0, diag, perm, diagShifted, shift);
#endif #endif
@ -914,7 +918,7 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
std::cout << "f(" << rightShifted << ") =" << fRight << " ; " << left << " " << shift << " " << right << "\n"; std::cout << "f(" << rightShifted << ") =" << fRight << " ; " << left << " " << shift << " " << right << "\n";
// assert((numext::isfinite)(fRight)); // assert((numext::isfinite)(fRight));
#endif #endif
#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
if(!(fLeft * fRight<0)) if(!(fLeft * fRight<0))
{ {
@ -948,7 +952,7 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
} }
muCur = (leftShifted + rightShifted) / Literal(2); muCur = (leftShifted + rightShifted) / Literal(2);
} }
else else
{ {
// We have a problem as shifting on the left or right give either a positive or negative value // We have a problem as shifting on the left or right give either a positive or negative value
// at the middle of [left,right]... // at the middle of [left,right]...
@ -959,7 +963,7 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
muCur = -muCur; muCur = -muCur;
} }
} }
singVals[k] = shift + muCur; singVals[k] = shift + muCur;
shifts[k] = shift; shifts[k] = shift;
mus[k] = muCur; mus[k] = muCur;
@ -1070,7 +1074,7 @@ void BDCSVD<MatrixType>::computeSingVecs
{ {
Index n = zhat.size(); Index n = zhat.size();
Index m = perm.size(); Index m = perm.size();
for (Index k = 0; k < n; ++k) for (Index k = 0; k < n; ++k)
{ {
if (zhat(k) == Literal(0)) if (zhat(k) == Literal(0))
@ -1088,7 +1092,7 @@ void BDCSVD<MatrixType>::computeSingVecs
} }
U(n,k) = Literal(0); U(n,k) = Literal(0);
U.col(k).normalize(); U.col(k).normalize();
if (m_compV) if (m_compV)
{ {
V.col(k).setZero(); V.col(k).setZero();
@ -1124,10 +1128,10 @@ void BDCSVD<MatrixType>::deflation43(Eigen::Index firstCol, Eigen::Index shift,
m_computed(start+i, start+i) = Literal(0); m_computed(start+i, start+i) = Literal(0);
return; return;
} }
m_computed(start,start) = r; m_computed(start,start) = r;
m_computed(start+i, start) = Literal(0); m_computed(start+i, start) = Literal(0);
m_computed(start+i, start+i) = Literal(0); m_computed(start+i, start+i) = Literal(0);
JacobiRotation<RealScalar> J(c/r,-s/r); JacobiRotation<RealScalar> J(c/r,-s/r);
if (m_compU) m_naiveU.middleRows(firstCol, size+1).applyOnTheRight(firstCol, firstCol+i, J); if (m_compU) m_naiveU.middleRows(firstCol, size+1).applyOnTheRight(firstCol, firstCol+i, J);
else m_naiveU.applyOnTheRight(firstCol, firstCol+i, J); else m_naiveU.applyOnTheRight(firstCol, firstCol+i, J);
@ -1184,29 +1188,29 @@ void BDCSVD<MatrixType>::deflation(Eigen::Index firstCol, Eigen::Index lastCol,
using std::sqrt; using std::sqrt;
using std::abs; using std::abs;
const Index length = lastCol + 1 - firstCol; const Index length = lastCol + 1 - firstCol;
Block<MatrixXr,Dynamic,1> col0(m_computed, firstCol+shift, firstCol+shift, length, 1); Block<MatrixXr,Dynamic,1> col0(m_computed, firstCol+shift, firstCol+shift, length, 1);
Diagonal<MatrixXr> fulldiag(m_computed); Diagonal<MatrixXr> fulldiag(m_computed);
VectorBlock<Diagonal<MatrixXr>,Dynamic> diag(fulldiag, firstCol+shift, length); VectorBlock<Diagonal<MatrixXr>,Dynamic> diag(fulldiag, firstCol+shift, length);
const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)(); const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)();
RealScalar maxDiag = diag.tail((std::max)(Index(1),length-1)).cwiseAbs().maxCoeff(); RealScalar maxDiag = diag.tail((std::max)(Index(1),length-1)).cwiseAbs().maxCoeff();
RealScalar epsilon_strict = numext::maxi<RealScalar>(considerZero,NumTraits<RealScalar>::epsilon() * maxDiag); RealScalar epsilon_strict = numext::maxi<RealScalar>(considerZero,NumTraits<RealScalar>::epsilon() * maxDiag);
RealScalar epsilon_coarse = Literal(8) * NumTraits<RealScalar>::epsilon() * numext::maxi<RealScalar>(col0.cwiseAbs().maxCoeff(), maxDiag); RealScalar epsilon_coarse = Literal(8) * NumTraits<RealScalar>::epsilon() * numext::maxi<RealScalar>(col0.cwiseAbs().maxCoeff(), maxDiag);
#ifdef EIGEN_BDCSVD_SANITY_CHECKS #ifdef EIGEN_BDCSVD_SANITY_CHECKS
assert(m_naiveU.allFinite()); assert(m_naiveU.allFinite());
assert(m_naiveV.allFinite()); assert(m_naiveV.allFinite());
assert(m_computed.allFinite()); assert(m_computed.allFinite());
#endif #endif
#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
std::cout << "\ndeflate:" << diag.head(k+1).transpose() << " | " << diag.segment(k+1,length-k-1).transpose() << "\n"; std::cout << "\ndeflate:" << diag.head(k+1).transpose() << " | " << diag.segment(k+1,length-k-1).transpose() << "\n";
#endif #endif
//condition 4.1 //condition 4.1
if (diag(0) < epsilon_coarse) if (diag(0) < epsilon_coarse)
{ {
#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
std::cout << "deflation 4.1, because " << diag(0) << " < " << epsilon_coarse << "\n"; std::cout << "deflation 4.1, because " << diag(0) << " < " << epsilon_coarse << "\n";
#endif #endif
@ -1243,22 +1247,22 @@ void BDCSVD<MatrixType>::deflation(Eigen::Index firstCol, Eigen::Index lastCol,
std::cout << " : " << col0.transpose() << "\n\n"; std::cout << " : " << col0.transpose() << "\n\n";
#endif #endif
{ {
// Check for total deflation // 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 // 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()<considerZero).all(); const bool total_deflation = (col0.tail(length-1).array().abs()<considerZero).all();
// Sort the diagonal entries, since diag(1:k-1) and diag(k:length) are already sorted, let's do a sorted merge. // 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. // First, compute the respective permutation.
Index *permutation = m_workspaceI.data(); Index *permutation = m_workspaceI.data();
{ {
permutation[0] = 0; permutation[0] = 0;
Index p = 1; Index p = 1;
// Move deflated diagonal entries at the end. // Move deflated diagonal entries at the end.
for(Index i=1; i<length; ++i) for(Index i=1; i<length; ++i)
if(abs(diag(i))<considerZero) if(abs(diag(i))<considerZero)
permutation[p++] = i; permutation[p++] = i;
Index i=1, j=k+1; Index i=1, j=k+1;
for( ; p < length; ++p) for( ; p < length; ++p)
{ {
@ -1268,7 +1272,7 @@ void BDCSVD<MatrixType>::deflation(Eigen::Index firstCol, Eigen::Index lastCol,
else permutation[p] = i++; else permutation[p] = i++;
} }
} }
// If we have a total deflation, then we have to insert diag(0) at the right place // If we have a total deflation, then we have to insert diag(0) at the right place
if(total_deflation) if(total_deflation)
{ {
@ -1284,22 +1288,22 @@ void BDCSVD<MatrixType>::deflation(Eigen::Index firstCol, Eigen::Index lastCol,
} }
} }
} }
// Current index of each col, and current column of each index // Current index of each col, and current column of each index
Index *realInd = m_workspaceI.data()+length; Index *realInd = m_workspaceI.data()+length;
Index *realCol = m_workspaceI.data()+2*length; Index *realCol = m_workspaceI.data()+2*length;
for(int pos = 0; pos< length; pos++) for(int pos = 0; pos< length; pos++)
{ {
realCol[pos] = pos; realCol[pos] = pos;
realInd[pos] = pos; realInd[pos] = pos;
} }
for(Index i = total_deflation?0:1; i < length; i++) for(Index i = total_deflation?0:1; i < length; i++)
{ {
const Index pi = permutation[length - (total_deflation ? i+1 : i)]; const Index pi = permutation[length - (total_deflation ? i+1 : i)];
const Index J = realCol[pi]; const Index J = realCol[pi];
using std::swap; using std::swap;
// swap diagonal and first column entries: // swap diagonal and first column entries:
swap(diag(i), diag(J)); swap(diag(i), diag(J));
@ -1322,7 +1326,7 @@ void BDCSVD<MatrixType>::deflation(Eigen::Index firstCol, Eigen::Index lastCol,
std::cout << "sorted: " << diag.transpose().format(bdcsvdfmt) << "\n"; std::cout << "sorted: " << diag.transpose().format(bdcsvdfmt) << "\n";
std::cout << " : " << col0.transpose() << "\n\n"; std::cout << " : " << col0.transpose() << "\n\n";
#endif #endif
//condition 4.4 //condition 4.4
{ {
Index i = length-1; Index i = length-1;
@ -1337,12 +1341,12 @@ void BDCSVD<MatrixType>::deflation(Eigen::Index firstCol, Eigen::Index lastCol,
deflation44(firstCol, firstCol + shift, firstRowW, firstColW, i-1, i, length); deflation44(firstCol, firstCol + shift, firstRowW, firstColW, i-1, i, length);
} }
} }
#ifdef EIGEN_BDCSVD_SANITY_CHECKS #ifdef EIGEN_BDCSVD_SANITY_CHECKS
for(Index j=2;j<length;++j) for(Index j=2;j<length;++j)
assert(diag(j-1)<=diag(j) || abs(diag(j))<considerZero); assert(diag(j-1)<=diag(j) || abs(diag(j))<considerZero);
#endif #endif
#ifdef EIGEN_BDCSVD_SANITY_CHECKS #ifdef EIGEN_BDCSVD_SANITY_CHECKS
assert(m_naiveU.allFinite()); assert(m_naiveU.allFinite());
assert(m_naiveV.allFinite()); assert(m_naiveV.allFinite());

51
test/bdcsvd.cpp
View File

@ -54,27 +54,53 @@ void bdcsvd_method()
VERIFY_IS_APPROX(m.bdcSvd(ComputeFullU|ComputeFullV).adjoint().solve(m), m); VERIFY_IS_APPROX(m.bdcSvd(ComputeFullU|ComputeFullV).adjoint().solve(m), m);
} }
// compare the Singular values returned with Jacobi and Bdc // Compare the Singular values returned with Jacobi and Bdc.
template<typename MatrixType> template<typename MatrixType>
void compare_bdc_jacobi(const MatrixType& a = MatrixType(), unsigned int computationOptions = 0) void compare_bdc_jacobi(const MatrixType& a = MatrixType(), unsigned int computationOptions = 0, int algoswap = 16, bool random = true)
{ {
MatrixType m = MatrixType::Random(a.rows(), a.cols()); MatrixType m = random ? MatrixType::Random(a.rows(), a.cols()) : a;
BDCSVD<MatrixType> bdc_svd(m);
BDCSVD<MatrixType> bdc_svd(m.rows(), m.cols(), computationOptions);
bdc_svd.setSwitchSize(algoswap);
bdc_svd.compute(m);
JacobiSVD<MatrixType> jacobi_svd(m); JacobiSVD<MatrixType> jacobi_svd(m);
VERIFY_IS_APPROX(bdc_svd.singularValues(), jacobi_svd.singularValues()); VERIFY_IS_APPROX(bdc_svd.singularValues(), jacobi_svd.singularValues());
if(computationOptions & ComputeFullU) VERIFY_IS_APPROX(bdc_svd.matrixU(), jacobi_svd.matrixU()); if(computationOptions & ComputeFullU) VERIFY_IS_APPROX(bdc_svd.matrixU(), jacobi_svd.matrixU());
if(computationOptions & ComputeThinU) VERIFY_IS_APPROX(bdc_svd.matrixU(), jacobi_svd.matrixU()); if(computationOptions & ComputeThinU) VERIFY_IS_APPROX(bdc_svd.matrixU(), jacobi_svd.matrixU());
if(computationOptions & ComputeFullV) VERIFY_IS_APPROX(bdc_svd.matrixV(), jacobi_svd.matrixV()); if(computationOptions & ComputeFullV) VERIFY_IS_APPROX(bdc_svd.matrixV(), jacobi_svd.matrixV());
if(computationOptions & ComputeThinV) VERIFY_IS_APPROX(bdc_svd.matrixV(), jacobi_svd.matrixV()); if(computationOptions & ComputeThinV) VERIFY_IS_APPROX(bdc_svd.matrixV(), jacobi_svd.matrixV());
} }
// Verifies total deflation is **not** triggered.
void compare_bdc_jacobi_instance(bool structure_as_m, int algoswap = 16)
{
MatrixXd m(4, 3);
if (structure_as_m) {
// The first 3 rows are the reduced form of Matrix 1 as shown below, and it
// has nonzero elements in the first column and diagonals only.
m << 1.056293, 0, 0,
-0.336468, 0.907359, 0,
-1.566245, 0, 0.149150,
-0.1, 0, 0;
} else {
// Matrix 1.
m << 0.882336, 18.3914, -26.7921,
-5.58135, 17.1931, -24.0892,
-20.794, 8.68496, -4.83103,
-8.4981, -10.5451, 23.9072;
}
compare_bdc_jacobi(m, 0, algoswap, false);
}
EIGEN_DECLARE_TEST(bdcsvd) EIGEN_DECLARE_TEST(bdcsvd)
{ {
CALL_SUBTEST_3(( svd_verify_assert<BDCSVD<Matrix3f> >(Matrix3f()) )); CALL_SUBTEST_3(( svd_verify_assert<BDCSVD<Matrix3f> >(Matrix3f()) ));
CALL_SUBTEST_4(( svd_verify_assert<BDCSVD<Matrix4d> >(Matrix4d()) )); CALL_SUBTEST_4(( svd_verify_assert<BDCSVD<Matrix4d> >(Matrix4d()) ));
CALL_SUBTEST_7(( svd_verify_assert<BDCSVD<MatrixXf> >(MatrixXf(10,12)) )); CALL_SUBTEST_7(( svd_verify_assert<BDCSVD<MatrixXf> >(MatrixXf(10,12)) ));
CALL_SUBTEST_8(( svd_verify_assert<BDCSVD<MatrixXcd> >(MatrixXcd(7,5)) )); CALL_SUBTEST_8(( svd_verify_assert<BDCSVD<MatrixXcd> >(MatrixXcd(7,5)) ));
CALL_SUBTEST_101(( svd_all_trivial_2x2(bdcsvd<Matrix2cd>) )); CALL_SUBTEST_101(( svd_all_trivial_2x2(bdcsvd<Matrix2cd>) ));
CALL_SUBTEST_102(( svd_all_trivial_2x2(bdcsvd<Matrix2d>) )); CALL_SUBTEST_102(( svd_all_trivial_2x2(bdcsvd<Matrix2d>) ));
@ -85,10 +111,10 @@ EIGEN_DECLARE_TEST(bdcsvd)
int r = internal::random<int>(1, EIGEN_TEST_MAX_SIZE/2), int r = internal::random<int>(1, EIGEN_TEST_MAX_SIZE/2),
c = internal::random<int>(1, EIGEN_TEST_MAX_SIZE/2); c = internal::random<int>(1, EIGEN_TEST_MAX_SIZE/2);
TEST_SET_BUT_UNUSED_VARIABLE(r) TEST_SET_BUT_UNUSED_VARIABLE(r)
TEST_SET_BUT_UNUSED_VARIABLE(c) TEST_SET_BUT_UNUSED_VARIABLE(c)
CALL_SUBTEST_6(( bdcsvd(Matrix<double,Dynamic,2>(r,2)) )); CALL_SUBTEST_6(( bdcsvd(Matrix<double,Dynamic,2>(r,2)) ));
CALL_SUBTEST_7(( bdcsvd(MatrixXf(r,c)) )); CALL_SUBTEST_7(( bdcsvd(MatrixXf(r,c)) ));
CALL_SUBTEST_7(( compare_bdc_jacobi(MatrixXf(r,c)) )); CALL_SUBTEST_7(( compare_bdc_jacobi(MatrixXf(r,c)) ));
@ -114,5 +140,12 @@ EIGEN_DECLARE_TEST(bdcsvd)
// CALL_SUBTEST_9( svd_preallocate<void>() ); // CALL_SUBTEST_9( svd_preallocate<void>() );
CALL_SUBTEST_2( svd_underoverflow<void>() ); CALL_SUBTEST_2( svd_underoverflow<void>() );
}
// Without total deflation issues.
CALL_SUBTEST_11(( compare_bdc_jacobi_instance(true) ));
CALL_SUBTEST_12(( compare_bdc_jacobi_instance(false) ));
// With total deflation issues before, when it shouldn't be triggered.
CALL_SUBTEST_13(( compare_bdc_jacobi_instance(true, 3) ));
CALL_SUBTEST_14(( compare_bdc_jacobi_instance(false, 3) ));
}