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

Fix BDCSVD's total deflation in branch 3.4, similar to that of master in MR 707.

Browse Source 
(cherry picked from commit 4d045eba53f9a32d052eb942448ba62def066529)
This commit is contained in:
Xinle Liu 2021-11-03 10:56:14 -07:00
parent 6b6ba41269
commit 9c193db5c7
2 changed files with 47 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -27,6 +27,10 @@
#define eigen_internal_assert(X) assert(X);
#endif
#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
#include <iostream>
#endif
namespace Eigen {
#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
@ -172,7 +176,7 @@ public:
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;
}
@ -404,7 +408,7 @@ void BDCSVD<MatrixType>::structured_update(Block<MatrixXr,Dynamic,Dynamic> A, co
//@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.
//@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
// 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>
@ -899,7 +903,7 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
RealScalar fLeft = secularEq(leftShifted, col0, diag, perm, diagShifted, shift);
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);
#endif
@ -1242,8 +1246,8 @@ void BDCSVD<MatrixType>::deflation(Eigen::Index firstCol, Eigen::Index lastCol,
#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()<considerZero).all();
// If we have a total deflation, then we have to consider col0(0)==diag(0) as a singular value during sorting.
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.
// First, compute the respective permutation.

42
test/bdcsvd.cpp
View File

@ -54,20 +54,46 @@ void bdcsvd_method()
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>
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());
BDCSVD<MatrixType> bdc_svd(m);
MatrixType m = random ? MatrixType::Random(a.rows(), a.cols()) : a;
BDCSVD<MatrixType> bdc_svd(m.rows(), m.cols(), computationOptions);
bdc_svd.setSwitchSize(algoswap);
bdc_svd.compute(m);
JacobiSVD<MatrixType> jacobi_svd(m);
VERIFY_IS_APPROX(bdc_svd.singularValues(), jacobi_svd.singularValues());
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 & ComputeFullV) 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)
{
CALL_SUBTEST_3(( svd_verify_assert<BDCSVD<Matrix3f> >(Matrix3f()) ));
@ -114,5 +140,13 @@ EIGEN_DECLARE_TEST(bdcsvd)
// CALL_SUBTEST_9( svd_preallocate<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) ));
}