From e093b43b2c40f00495937c3134bf55ba29676993 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Mon, 6 Jul 2009 17:12:10 +0200 Subject: [PATCH 1/3] * rename QR to HouseholderQR because really that impacts the API, not just the impl. * rename qr() to householderQr(), for same reason. * clarify that it's non-pivoting, non-rank-revealing, so remove all the rank API, make solve() be void instead of bool, update the docs/test, etc. * fix warning in SVD --- Eigen/src/Core/MatrixBase.h | 2 +- Eigen/src/Core/util/ForwardDeclarations.h | 2 +- Eigen/src/QR/QR.h | 161 ++++-------------- Eigen/src/SVD/SVD.h | 1 - doc/eigendoxy_header.html.in | 2 +- .../{QR_solve.cpp => HouseholderQR_solve.cpp} | 11 +- test/geo_hyperplane.cpp | 2 +- test/main.h | 4 +- test/qr.cpp | 62 +------ 9 files changed, 43 insertions(+), 204 deletions(-) rename doc/snippets/{QR_solve.cpp => HouseholderQR_solve.cpp} (52%) diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 65ab02d62..941539214 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -653,7 +653,7 @@ template class MatrixBase /////////// QR module /////////// - const QR qr() const; + const HouseholderQR householderQr() const; EigenvaluesReturnType eigenvalues() const; RealScalar operatorNorm() const; diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index b457268af..a2105604a 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -112,7 +112,7 @@ template class Reverse; template class LU; template class PartialLU; -template class QR; +template class HouseholderQR; template class SVD; template class LLT; template class LDLT; diff --git a/Eigen/src/QR/QR.h b/Eigen/src/QR/QR.h index 90751dd42..d6d3d2081 100644 --- a/Eigen/src/QR/QR.h +++ b/Eigen/src/QR/QR.h @@ -22,24 +22,26 @@ // License and a copy of the GNU General Public License along with // Eigen. If not, see . -#ifndef EIGEN_QR_H -#define EIGEN_QR_H +#ifndef EIGEN_HouseholderQR_H +#define EIGEN_HouseholderQR_H -/** \ingroup QR_Module +/** \ingroup HouseholderQR_Module * \nonstableyet * - * \class QR + * \class HouseholderQR * - * \brief QR decomposition of a matrix + * \brief Householder QR decomposition of a matrix * * \param MatrixType the type of the matrix of which we are computing the QR decomposition * * This class performs a QR decomposition using Householder transformations. The result is * stored in a compact way compatible with LAPACK. * + * Note that no pivoting is performed. This is \b not a rank-revealing decomposition. + * * \sa MatrixBase::qr() */ -template class QR +template class HouseholderQR { public: @@ -53,88 +55,23 @@ template class QR * \brief Default Constructor. * * The default constructor is useful in cases in which the user intends to - * perform decompositions via QR::compute(const MatrixType&). + * perform decompositions via HouseholderQR::compute(const MatrixType&). */ - QR() : m_qr(), m_hCoeffs(), m_isInitialized(false) {} + HouseholderQR() : m_qr(), m_hCoeffs(), m_isInitialized(false) {} - QR(const MatrixType& matrix) + HouseholderQR(const MatrixType& matrix) : m_qr(matrix.rows(), matrix.cols()), m_hCoeffs(matrix.cols()), m_isInitialized(false) { compute(matrix); } - - /** \deprecated use isInjective() - * \returns whether or not the matrix is of full rank - * - * \note Since the rank is computed only once, i.e. the first time it is needed, this - * method almost does not perform any further computation. - */ - EIGEN_DEPRECATED bool isFullRank() const - { - ei_assert(m_isInitialized && "QR is not initialized."); - return rank() == m_qr.cols(); - } - - /** \returns the rank of the matrix of which *this is the QR decomposition. - * - * \note Since the rank is computed only once, i.e. the first time it is needed, this - * method almost does not perform any further computation. - */ - int rank() const; - - /** \returns the dimension of the kernel of the matrix of which *this is the QR decomposition. - * - * \note Since the rank is computed only once, i.e. the first time it is needed, this - * method almost does not perform any further computation. - */ - inline int dimensionOfKernel() const - { - ei_assert(m_isInitialized && "QR is not initialized."); - return m_qr.cols() - rank(); - } - - /** \returns true if the matrix of which *this is the QR decomposition represents an injective - * linear map, i.e. has trivial kernel; false otherwise. - * - * \note Since the rank is computed only once, i.e. the first time it is needed, this - * method almost does not perform any further computation. - */ - inline bool isInjective() const - { - ei_assert(m_isInitialized && "QR is not initialized."); - return rank() == m_qr.cols(); - } - - /** \returns true if the matrix of which *this is the QR decomposition represents a surjective - * linear map; false otherwise. - * - * \note Since the rank is computed only once, i.e. the first time it is needed, this - * method almost does not perform any further computation. - */ - inline bool isSurjective() const - { - ei_assert(m_isInitialized && "QR is not initialized."); - return rank() == m_qr.rows(); - } - - /** \returns true if the matrix of which *this is the QR decomposition is invertible. - * - * \note Since the rank is computed only once, i.e. the first time it is needed, this - * method almost does not perform any further computation. - */ - inline bool isInvertible() const - { - ei_assert(m_isInitialized && "QR is not initialized."); - return isInjective() && isSurjective(); - } - + /** \returns a read-only expression of the matrix R of the actual the QR decomposition */ const Part, UpperTriangular> matrixR(void) const { - ei_assert(m_isInitialized && "QR is not initialized."); + ei_assert(m_isInitialized && "HouseholderQR is not initialized."); int cols = m_qr.cols(); return MatrixRBlockType(m_qr, 0, 0, cols, cols).nestByValue().template part(); } @@ -148,22 +85,14 @@ template class QR * Resized if necessary, so that result->rows()==A.cols() and result->cols()==b.cols(). * If no solution exists, *result is left with undefined coefficients. * - * \returns true if any solution exists, false if no solution exists. - * - * \note If there exist more than one solution, this method will arbitrarily choose one. - * If you need a complete analysis of the space of solutions, take the one solution obtained - * by this method and add to it elements of the kernel, as determined by kernel(). - * * \note The case where b is a matrix is not yet implemented. Also, this * code is space inefficient. * - * Example: \include QR_solve.cpp - * Output: \verbinclude QR_solve.out - * - * \sa MatrixBase::solveTriangular(), kernel(), computeKernel(), inverse(), computeInverse() + * Example: \include HouseholderQR_solve.cpp + * Output: \verbinclude HouseholderQR_solve.out */ template - bool solve(const MatrixBase& b, ResultType *result) const; + void solve(const MatrixBase& b, ResultType *result) const; MatrixType matrixQ(void) const; @@ -172,34 +101,14 @@ template class QR protected: MatrixType m_qr; VectorType m_hCoeffs; - mutable int m_rank; - mutable bool m_rankIsUptodate; bool m_isInitialized; }; -/** \returns the rank of the matrix of which *this is the QR decomposition. */ -template -int QR::rank() const -{ - ei_assert(m_isInitialized && "QR is not initialized."); - if (!m_rankIsUptodate) - { - RealScalar maxCoeff = m_qr.diagonal().cwise().abs().maxCoeff(); - int n = m_qr.cols(); - m_rank = 0; - while(m_rank -void QR::compute(const MatrixType& matrix) +void HouseholderQR::compute(const MatrixType& matrix) { - m_rankIsUptodate = false; m_qr = matrix; m_hCoeffs.resize(matrix.cols()); @@ -262,12 +171,12 @@ void QR::compute(const MatrixType& matrix) template template -bool QR::solve( +void HouseholderQR::solve( const MatrixBase& b, ResultType *result ) const { - ei_assert(m_isInitialized && "QR is not initialized."); + ei_assert(m_isInitialized && "HouseholderQR is not initialized."); const int rows = m_qr.rows(); ei_assert(b.rows() == rows); result->resize(rows, b.cols()); @@ -276,27 +185,17 @@ bool QR::solve( // Q^T without explicitly forming matrixQ(). Investigate. *result = matrixQ().transpose()*b; - if(!isSurjective()) - { - // is result is in the image of R ? - RealScalar biggest_in_res = result->corner(TopLeft, m_rank, result->cols()).cwise().abs().maxCoeff(); - for(int col = 0; col < result->cols(); ++col) - for(int row = m_rank; row < result->rows(); ++row) - if(!ei_isMuchSmallerThan(result->coeff(row,col), biggest_in_res)) - return false; - } - m_qr.corner(TopLeft, m_rank, m_rank) + const int rank = std::min(result->rows(), result->cols()); + m_qr.corner(TopLeft, rank, rank) .template marked() - .solveTriangularInPlace(result->corner(TopLeft, m_rank, result->cols())); - - return true; + .solveTriangularInPlace(result->corner(TopLeft, rank, result->cols())); } /** \returns the matrix Q */ template -MatrixType QR::matrixQ() const +MatrixType HouseholderQR::matrixQ() const { - ei_assert(m_isInitialized && "QR is not initialized."); + ei_assert(m_isInitialized && "HouseholderQR is not initialized."); // compute the product Q_0 Q_1 ... Q_n-1, // where Q_k is the k-th Householder transformation I - h_k v_k v_k' // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...] @@ -319,16 +218,16 @@ MatrixType QR::matrixQ() const #endif // EIGEN_HIDE_HEAVY_CODE -/** \return the QR decomposition of \c *this. +/** \return the Householder QR decomposition of \c *this. * - * \sa class QR + * \sa class HouseholderQR */ template -const QR::PlainMatrixType> -MatrixBase::qr() const +const HouseholderQR::PlainMatrixType> +MatrixBase::householderQr() const { - return QR(eval()); + return HouseholderQR(eval()); } -#endif // EIGEN_QR_H +#endif // EIGEN_HouseholderQR_H diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h index 71f6763a8..97f82c78f 100644 --- a/Eigen/src/SVD/SVD.h +++ b/Eigen/src/SVD/SVD.h @@ -145,7 +145,6 @@ void SVD::compute(const MatrixType& matrix) { const int m = matrix.rows(); const int n = matrix.cols(); - const int nu = std::min(m,n); m_matU.resize(m, m); m_matU.setZero(); diff --git a/doc/eigendoxy_header.html.in b/doc/eigendoxy_header.html.in index 97cafd41d..572c47158 100644 --- a/doc/eigendoxy_header.html.in +++ b/doc/eigendoxy_header.html.in @@ -6,5 +6,5 @@ \ No newline at end of file diff --git a/doc/snippets/QR_solve.cpp b/doc/snippets/HouseholderQR_solve.cpp similarity index 52% rename from doc/snippets/QR_solve.cpp rename to doc/snippets/HouseholderQR_solve.cpp index 7e629f851..aa9404951 100644 --- a/doc/snippets/QR_solve.cpp +++ b/doc/snippets/HouseholderQR_solve.cpp @@ -4,11 +4,6 @@ Matrix3f y = Matrix3f::Random(); cout << "Here is the matrix m:" << endl << m << endl; cout << "Here is the matrix y:" << endl << y << endl; Matrix3f x; -if(m.qr().solve(y, &x)) -{ - assert(y.isApprox(m*x)); - cout << "Here is a solution x to the equation mx=y:" << endl << x << endl; -} -else - cout << "The equation mx=y does not have any solution." << endl; - +m.householderQr().solve(y, &x)); +assert(y.isApprox(m*x)); +cout << "Here is a solution x to the equation mx=y:" << endl << x << endl; diff --git a/test/geo_hyperplane.cpp b/test/geo_hyperplane.cpp index c5d2715db..f8281a16b 100644 --- a/test/geo_hyperplane.cpp +++ b/test/geo_hyperplane.cpp @@ -64,7 +64,7 @@ template void hyperplane(const HyperplaneType& _plane) // transform if (!NumTraits::IsComplex) { - MatrixType rot = MatrixType::Random(dim,dim).qr().matrixQ(); + MatrixType rot = MatrixType::Random(dim,dim).householderQr().matrixQ(); DiagonalMatrix scaling(VectorType::Random()); Translation translation(VectorType::Random()); diff --git a/test/main.h b/test/main.h index a6a780603..53c8245c6 100644 --- a/test/main.h +++ b/test/main.h @@ -241,8 +241,8 @@ void createRandomMatrixOfRank(int desired_rank, int rows, int cols, Eigen::Matri const int diag_size = std::min(d.rows(),d.cols()); d.diagonal().segment(desired_rank, diag_size-desired_rank) = VectorType::Zero(diag_size-desired_rank); - QR qra(a); - QR qrb(b); + HouseholderQR qra(a); + HouseholderQR qrb(b); m = (qra.matrixQ() * d * qrb.matrixQ()).lazy(); } diff --git a/test/qr.cpp b/test/qr.cpp index 6e96f1e97..88a447c4b 100644 --- a/test/qr.cpp +++ b/test/qr.cpp @@ -36,7 +36,7 @@ template void qr(const MatrixType& m) typedef Matrix VectorType; MatrixType a = MatrixType::Random(rows,cols); - QR qrOfA(a); + HouseholderQR qrOfA(a); VERIFY_IS_APPROX(a, qrOfA.matrixQ() * qrOfA.matrixR()); VERIFY_IS_NOT_APPROX(a+MatrixType::Identity(rows, cols), qrOfA.matrixQ() * qrOfA.matrixR()); @@ -55,40 +55,6 @@ template void qr(const MatrixType& m) VERIFY_IS_APPROX(b, hess.matrixQ() * hess.matrixH() * hess.matrixQ().adjoint()); } -template void qr_non_invertible() -{ - /* this test covers the following files: QR.h */ - int rows = ei_random(20,200), cols = ei_random(20,rows), cols2 = ei_random(20,rows); - int rank = ei_random(1, std::min(rows, cols)-1); - - MatrixType m1(rows, cols), m2(cols, cols2), m3(rows, cols2), k(1,1); - createRandomMatrixOfRank(rank, rows, cols, m1); - - QR lu(m1); -// typename LU::KernelResultType m1kernel = lu.kernel(); -// typename LU::ImageResultType m1image = lu.image(); - std::cerr << rows << "x" << cols << " " << rank << " " << lu.rank() << "\n"; - if (rank != lu.rank()) - std::cerr << lu.matrixR().diagonal().transpose() << "\n"; - VERIFY(rank == lu.rank()); - VERIFY(cols - lu.rank() == lu.dimensionOfKernel()); - VERIFY(!lu.isInjective()); - VERIFY(!lu.isInvertible()); - VERIFY(lu.isSurjective() == (lu.rank() == rows)); -// VERIFY((m1 * m1kernel).isMuchSmallerThan(m1)); -// VERIFY(m1image.lu().rank() == rank); -// MatrixType sidebyside(m1.rows(), m1.cols() + m1image.cols()); -// sidebyside << m1, m1image; -// VERIFY(sidebyside.lu().rank() == rank); - m2 = MatrixType::Random(cols,cols2); - m3 = m1*m2; - m2 = MatrixType::Random(cols,cols2); - lu.solve(m3, &m2); - VERIFY_IS_APPROX(m3, m1*m2); - m3 = MatrixType::Random(rows,cols2); - VERIFY(!lu.solve(m3, &m2)); -} - template void qr_invertible() { /* this test covers the following files: QR.h */ @@ -105,33 +71,18 @@ template void qr_invertible() m1 += a * a.adjoint(); } - QR lu(m1); - VERIFY(0 == lu.dimensionOfKernel()); - VERIFY(size == lu.rank()); - VERIFY(lu.isInjective()); - VERIFY(lu.isSurjective()); - VERIFY(lu.isInvertible()); -// VERIFY(lu.image().lu().isInvertible()); + HouseholderQR qr(m1); m3 = MatrixType::Random(size,size); - lu.solve(m3, &m2); + qr.solve(m3, &m2); //std::cerr << m3 - m1*m2 << "\n\n"; VERIFY_IS_APPROX(m3, m1*m2); -// VERIFY_IS_APPROX(m2, lu.inverse()*m3); - m3 = MatrixType::Random(size,size); - VERIFY(lu.solve(m3, &m2)); } template void qr_verify_assert() { MatrixType tmp; - QR qr; - VERIFY_RAISES_ASSERT(qr.isFullRank()) - VERIFY_RAISES_ASSERT(qr.rank()) - VERIFY_RAISES_ASSERT(qr.dimensionOfKernel()) - VERIFY_RAISES_ASSERT(qr.isInjective()) - VERIFY_RAISES_ASSERT(qr.isSurjective()) - VERIFY_RAISES_ASSERT(qr.isInvertible()) + HouseholderQR qr; VERIFY_RAISES_ASSERT(qr.matrixR()) VERIFY_RAISES_ASSERT(qr.solve(tmp,&tmp)) VERIFY_RAISES_ASSERT(qr.matrixQ()) @@ -149,11 +100,6 @@ void test_qr() } for(int i = 0; i < g_repeat; i++) { - CALL_SUBTEST( qr_non_invertible() ); - CALL_SUBTEST( qr_non_invertible() ); - // TODO fix issue with complex -// CALL_SUBTEST( qr_non_invertible() ); -// CALL_SUBTEST( qr_non_invertible() ); CALL_SUBTEST( qr_invertible() ); CALL_SUBTEST( qr_invertible() ); // TODO fix issue with complex From e057beee4ee682f3128d36c74f3aa57fc9ce148e Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Mon, 6 Jul 2009 17:20:07 +0200 Subject: [PATCH 2/3] fix some search-and-replace damage --- Eigen/src/QR/QR.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Eigen/src/QR/QR.h b/Eigen/src/QR/QR.h index d6d3d2081..a83ec5ec8 100644 --- a/Eigen/src/QR/QR.h +++ b/Eigen/src/QR/QR.h @@ -22,10 +22,10 @@ // License and a copy of the GNU General Public License along with // Eigen. If not, see . -#ifndef EIGEN_HouseholderQR_H -#define EIGEN_HouseholderQR_H +#ifndef EIGEN_QR_H +#define EIGEN_QR_H -/** \ingroup HouseholderQR_Module +/** \ingroup QR_Module * \nonstableyet * * \class HouseholderQR @@ -230,4 +230,4 @@ MatrixBase::householderQr() const } -#endif // EIGEN_HouseholderQR_H +#endif // EIGEN_QR_H From 889726bf7cf82b30e4a7140f467a175bab1dca2d Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Mon, 6 Jul 2009 17:24:11 +0200 Subject: [PATCH 3/3] add matrixQR() method exposing the storage. that's where the householder thing impacts the API. --- Eigen/src/QR/QR.h | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/Eigen/src/QR/QR.h b/Eigen/src/QR/QR.h index a83ec5ec8..c5e8f9097 100644 --- a/Eigen/src/QR/QR.h +++ b/Eigen/src/QR/QR.h @@ -95,6 +95,11 @@ template class HouseholderQR void solve(const MatrixBase& b, ResultType *result) const; MatrixType matrixQ(void) const; + + /** \returns a reference to the matrix where the Householder QR decomposition is stored + * in a LAPACK-compatible way. + */ + const MatrixType& matrixQR() const { return m_qr; } void compute(const MatrixType& matrix);