From d92df336ad21c7f8e0289f8ac3084b6313a17fe4 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Tue, 23 Feb 2010 15:40:24 -0500 Subject: [PATCH] Further LU test improvements. I'm not aware of any test failures anymore, not even with huge numbers of repetitions. Finally the createRandomMatrixOfRank() function is renamed to createRandomPIMatrixOfRank, where PI stands for 'partial isometry', that is, a matrix whose singular values are 0 or 1. --- Eigen/src/LU/FullPivLU.h | 7 +++-- test/inverse.cpp | 2 +- test/lu.cpp | 65 +++++++++++++++------------------------- test/main.h | 12 ++++++-- test/qr_colpivoting.cpp | 4 +-- test/qr_fullpivoting.cpp | 2 +- 6 files changed, 42 insertions(+), 50 deletions(-) diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h index ec551645b..0a305d52b 100644 --- a/Eigen/src/LU/FullPivLU.h +++ b/Eigen/src/LU/FullPivLU.h @@ -422,8 +422,11 @@ FullPivLU& FullPivLU::compute(const MatrixType& matrix) // when k==0, biggest_in_corner is the biggest coeff absolute value in the original matrix if(k == 0) cutoff = biggest_in_corner * NumTraits::epsilon(); - // if the pivot (hence the corner) is exactly zero, terminate to avoid generating nan/inf values - if(ei_abs(biggest_in_corner) < cutoff) + // if the pivot (hence the corner) is "zero", terminate to avoid generating nan/inf values. + // Notice that using an exact comparison (biggest_in_corner==0) here, as Golub-van Loan do in + // their pseudo-code, results in numerical instability! The cutoff here has been validated + // by running the unit test 'lu' with many repetitions. + if(biggest_in_corner < cutoff) { // before exiting, make sure to initialize the still uninitialized transpositions // in a sane state without destroying what we already have. diff --git a/test/inverse.cpp b/test/inverse.cpp index 3f6138e0c..1e567ad14 100644 --- a/test/inverse.cpp +++ b/test/inverse.cpp @@ -42,7 +42,7 @@ template void inverse(const MatrixType& m) m2(rows, cols), mzero = MatrixType::Zero(rows, cols), identity = MatrixType::Identity(rows, rows); - createRandomProjectionOfRank(rows,rows,rows,m1); + createRandomPIMatrixOfRank(rows,rows,rows,m1); m2 = m1.inverse(); VERIFY_IS_APPROX(m1, m2.inverse() ); diff --git a/test/lu.cpp b/test/lu.cpp index 02f6ec805..442202a33 100644 --- a/test/lu.cpp +++ b/test/lu.cpp @@ -28,9 +28,6 @@ using namespace std; template void lu_non_invertible() { - static int times_called = 0; - times_called++; - typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; /* this test covers the following files: @@ -55,11 +52,16 @@ template void lu_non_invertible() cols2 = cols = MatrixType::ColsAtCompileTime; } + enum { + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + ColsAtCompileTime = MatrixType::ColsAtCompileTime + }; typedef typename ei_kernel_retval_base >::ReturnType KernelMatrixType; typedef typename ei_image_retval_base >::ReturnType ImageMatrixType; - typedef Matrix DynamicMatrixType; - typedef Matrix + typedef Matrix CMatrixType; + typedef Matrix + RMatrixType; int rank = ei_random(1, std::min(rows, cols)-1); @@ -68,26 +70,21 @@ template void lu_non_invertible() MatrixType m1(rows, cols), m3(rows, cols2); CMatrixType m2(cols, cols2); - createRandomProjectionOfRank(rank, rows, cols, m1); + createRandomPIMatrixOfRank(rank, rows, cols, m1); FullPivLU lu; - // The special value 0.01 below works well in tests. Keep in mind that we're only computing the rank of projections. - // So it's not clear at all the epsilon should play any role there. + // The special value 0.01 below works well in tests. Keep in mind that we're only computing the rank + // of singular values are either 0 or 1. + // So it's not clear at all that the epsilon should play any role there. lu.setThreshold(RealScalar(0.01)); lu.compute(m1); - // FIXME need better way to construct trapezoid matrices. extend triangularView to support rectangular. - DynamicMatrixType u(rows,cols); - for(int i = 0; i < rows; i++) - for(int j = 0; j < cols; j++) - u(i,j) = i>j ? Scalar(0) : lu.matrixLU()(i,j); - DynamicMatrixType l(rows,rows); - for(int i = 0; i < rows; i++) - for(int j = 0; j < rows; j++) - l(i,j) = (i=cols) ? Scalar(0) - : i==j ? Scalar(1) - : lu.matrixLU()(i,j); + MatrixType u(rows,cols); + u = lu.matrixLU().template triangularView(); + RMatrixType l = RMatrixType::Identity(rows,rows); + l.block(0,0,rows,std::min(rows,cols)).template triangularView() + = lu.matrixLU().block(0,0,rows,std::min(rows,cols)); VERIFY_IS_APPROX(lu.permutationP() * m1 * lu.permutationQ(), l*u); @@ -101,20 +98,8 @@ template void lu_non_invertible() VERIFY(!lu.isSurjective()); VERIFY((m1 * m1kernel).isMuchSmallerThan(m1)); VERIFY(m1image.fullPivLu().rank() == rank); + VERIFY_IS_APPROX(m1 * m1.adjoint() * m1image, m1image); - // The following test is damn hard to get to succeed over a large number of repetitions. - // We're checking that the image is indeed the image, i.e. adding it as new columns doesn't increase the rank. - // Since we've already tested rank() above, the point here is not to test rank(), it is to test image(). - // Since image() is implemented in a very simple way that doesn't leave much room for choice, the occasional - // errors that we get here (one in 1e+4 repetitions roughly) are probably just a sign that it's a really - // hard test, so we just limit how many times it's run. - if(times_called < 100) - { - DynamicMatrixType sidebyside(m1.rows(), m1.cols() + m1image.cols()); - sidebyside << m1, m1image; - VERIFY(sidebyside.fullPivLu().rank() == rank); - } - m2 = CMatrixType::Random(cols,cols2); m3 = m1*m2; m2 = CMatrixType::Random(cols,cols2); @@ -128,20 +113,18 @@ template void lu_invertible() /* this test covers the following files: LU.h */ + typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; int size = ei_random(1,200); MatrixType m1(size, size), m2(size, size), m3(size, size); - m1 = MatrixType::Random(size,size); + FullPivLU lu; + lu.setThreshold(0.01); + do { + m1 = MatrixType::Random(size,size); + lu.compute(m1); + } while(!lu.isInvertible()); - if (ei_is_same_type::ret) - { - // let's build a matrix more stable to inverse - MatrixType a = MatrixType::Random(size,size*2); - m1 += a * a.adjoint(); - } - - FullPivLU lu(m1); VERIFY(0 == lu.dimensionOfKernel()); VERIFY(lu.kernel().cols() == 1); // the kernel() should consist of a single (zero) column vector VERIFY(size == lu.rank()); diff --git a/test/main.h b/test/main.h index 6d296b2e3..96324de33 100644 --- a/test/main.h +++ b/test/main.h @@ -148,7 +148,7 @@ namespace Eigen #define EIGEN_INTERNAL_DEBUGGING #define EIGEN_NICE_RANDOM -#include // required for createRandomProjectionOfRank +#include // required for createRandomPIMatrixOfRank #define VERIFY(a) do { if (!(a)) { \ @@ -342,8 +342,13 @@ inline bool test_isUnitary(const MatrixBase& m) return m.isUnitary(test_precision::Scalar>()); } +/** Creates a random Partial Isometry matrix of given rank. + * + * A partial isometry is a matrix all of whose singular values are either 0 or 1. + * This is very useful to test rank-revealing algorithms. + */ template -void createRandomProjectionOfRank(int desired_rank, int rows, int cols, MatrixType& m) +void createRandomPIMatrixOfRank(int desired_rank, int rows, int cols, MatrixType& m) { typedef typename ei_traits::Scalar Scalar; enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime }; @@ -360,7 +365,8 @@ void createRandomProjectionOfRank(int desired_rank, int rows, int cols, MatrixTy if(desired_rank == 1) { - m = VectorType::Random(rows) * VectorType::Random(cols).transpose(); + // here we normalize the vectors to get a partial isometry + m = VectorType::Random(rows).normalized() * VectorType::Random(cols).normalized().transpose(); return; } diff --git a/test/qr_colpivoting.cpp b/test/qr_colpivoting.cpp index abee32184..96cc66316 100644 --- a/test/qr_colpivoting.cpp +++ b/test/qr_colpivoting.cpp @@ -36,7 +36,7 @@ template void qr() typedef Matrix MatrixQType; typedef Matrix VectorType; MatrixType m1; - createRandomProjectionOfRank(rank,rows,cols,m1); + createRandomPIMatrixOfRank(rank,rows,cols,m1); ColPivHouseholderQR qr(m1); VERIFY_IS_APPROX(rank, qr.rank()); VERIFY(cols - qr.rank() == qr.dimensionOfKernel()); @@ -64,7 +64,7 @@ template void qr_fixedsize() typedef typename MatrixType::Scalar Scalar; int rank = ei_random(1, std::min(int(Rows), int(Cols))-1); Matrix m1; - createRandomProjectionOfRank(rank,Rows,Cols,m1); + createRandomPIMatrixOfRank(rank,Rows,Cols,m1); ColPivHouseholderQR > qr(m1); VERIFY_IS_APPROX(rank, qr.rank()); VERIFY(Cols - qr.rank() == qr.dimensionOfKernel()); diff --git a/test/qr_fullpivoting.cpp b/test/qr_fullpivoting.cpp index 60255f94c..7ad3af1fe 100644 --- a/test/qr_fullpivoting.cpp +++ b/test/qr_fullpivoting.cpp @@ -35,7 +35,7 @@ template void qr() typedef Matrix MatrixQType; typedef Matrix VectorType; MatrixType m1; - createRandomProjectionOfRank(rank,rows,cols,m1); + createRandomPIMatrixOfRank(rank,rows,cols,m1); FullPivHouseholderQR qr(m1); VERIFY_IS_APPROX(rank, qr.rank()); VERIFY(cols - qr.rank() == qr.dimensionOfKernel());