diff --git a/Eigen/src/QR/EigenSolver.h b/Eigen/src/QR/EigenSolver.h index e9da5d941..70f21cebc 100644 --- a/Eigen/src/QR/EigenSolver.h +++ b/Eigen/src/QR/EigenSolver.h @@ -1,5 +1,5 @@ // This file is part of Eigen, a lightweight C++ template library -// for linear algebra. Eigen itself is part of the KDE project. +// for linear algebra. // // Copyright (C) 2008 Gael Guennebaud // @@ -53,9 +53,18 @@ template class EigenSolver typedef Matrix RealVectorType; typedef Matrix RealVectorTypeX; + /** + * \brief Default Constructor. + * + * The default constructor is useful in cases in which the user intends to + * perform decompositions via EigenSolver::compute(const MatrixType&). + */ + EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false) {} + EigenSolver(const MatrixType& matrix) : m_eivec(matrix.rows(), matrix.cols()), - m_eivalues(matrix.cols()) + m_eivalues(matrix.cols()), + m_isInitialized(false) { compute(matrix); } @@ -94,12 +103,20 @@ template class EigenSolver * * \sa pseudoEigenvalueMatrix() */ - const MatrixType& pseudoEigenvectors() const { return m_eivec; } + const MatrixType& pseudoEigenvectors() const + { + ei_assert(m_isInitialized && "EigenSolver is not initialized."); + return m_eivec; + } MatrixType pseudoEigenvalueMatrix() const; /** \returns the eigenvalues as a column vector */ - EigenvalueType eigenvalues() const { return m_eivalues; } + EigenvalueType eigenvalues() const + { + ei_assert(m_isInitialized && "EigenSolver is not initialized."); + return m_eivalues; + } void compute(const MatrixType& matrix); @@ -111,6 +128,7 @@ template class EigenSolver protected: MatrixType m_eivec; EigenvalueType m_eivalues; + bool m_isInitialized; }; /** \returns the real block diagonal matrix D of the eigenvalues. @@ -120,6 +138,7 @@ template class EigenSolver template MatrixType EigenSolver::pseudoEigenvalueMatrix() const { + ei_assert(m_isInitialized && "EigenSolver is not initialized."); int n = m_eivec.cols(); MatrixType matD = MatrixType::Zero(n,n); for (int i=0; i::pseudoEigenvalueMatrix() const template typename EigenSolver::EigenvectorType EigenSolver::eigenvectors(void) const { + ei_assert(m_isInitialized && "EigenSolver is not initialized."); int n = m_eivec.cols(); EigenvectorType matV(n,n); for (int j=0; j::compute(const MatrixType& matrix) // Reduce Hessenberg to real Schur form. hqr2(matH); + + m_isInitialized = true; } // Nonsymmetric reduction to Hessenberg form. diff --git a/Eigen/src/QR/QR.h b/Eigen/src/QR/QR.h index 13d49a4a3..90751dd42 100644 --- a/Eigen/src/QR/QR.h +++ b/Eigen/src/QR/QR.h @@ -1,5 +1,5 @@ // This file is part of Eigen, a lightweight C++ template library -// for linear algebra. Eigen itself is part of the KDE project. +// for linear algebra. // // Copyright (C) 2008 Gael Guennebaud // @@ -49,51 +49,146 @@ template class QR typedef Matrix MatrixTypeR; typedef Matrix VectorType; + /** + * \brief Default Constructor. + * + * The default constructor is useful in cases in which the user intends to + * perform decompositions via QR::compute(const MatrixType&). + */ + QR() : m_qr(), m_hCoeffs(), m_isInitialized(false) {} + QR(const MatrixType& matrix) : m_qr(matrix.rows(), matrix.cols()), - m_hCoeffs(matrix.cols()) + m_hCoeffs(matrix.cols()), + m_isInitialized(false) { - _compute(matrix); + 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 whether or not the matrix is of full rank */ - bool isFullRank() const { return rank() == std::min(m_qr.rows(),m_qr.cols()); } + /** \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(); + } - int rank() const; - /** \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."); int cols = m_qr.cols(); return MatrixRBlockType(m_qr, 0, 0, cols, cols).nestByValue().template part(); } + /** This method finds a solution x to the equation Ax=b, where A is the matrix of which + * *this is the QR decomposition, if any exists. + * + * \param b the right-hand-side of the equation to solve. + * + * \param result a pointer to the vector/matrix in which to store the solution, if any exists. + * 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() + */ + template + bool solve(const MatrixBase& b, ResultType *result) const; + MatrixType matrixQ(void) const; - private: - - void _compute(const MatrixType& matrix); + void compute(const MatrixType& matrix); 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().maxCoeff(); - int n = std::min(m_qr.rows(),m_qr.cols()); - m_rank = n; - for (int i=0; i::rank() const #ifndef EIGEN_HIDE_HEAVY_CODE template -void QR::_compute(const MatrixType& matrix) -{ +void QR::compute(const MatrixType& matrix) +{ m_rankIsUptodate = false; m_qr = matrix; + m_hCoeffs.resize(matrix.cols()); + int rows = matrix.rows(); int cols = matrix.cols(); + RealScalar eps2 = precision()*precision(); for (int k = 0; k < cols; ++k) { @@ -132,7 +230,8 @@ void QR::_compute(const MatrixType& matrix) m_hCoeffs.coeffRef(k) = 0; } } - else if ( (!ei_isMuchSmallerThan(beta=m_qr.col(k).end(remainingSize-1).squaredNorm(),static_cast(1))) || ei_imag(v0)==0 ) + else if ((beta=m_qr.col(k).end(remainingSize-1).squaredNorm())>eps2) + // FIXME what about ei_imag(v0) ?? { // form k-th Householder vector beta = ei_sqrt(ei_abs2(v0)+beta); @@ -158,12 +257,46 @@ void QR::_compute(const MatrixType& matrix) m_hCoeffs.coeffRef(k) = 0; } } + m_isInitialized = true; +} + +template +template +bool QR::solve( + const MatrixBase& b, + ResultType *result +) const +{ + ei_assert(m_isInitialized && "QR is not initialized."); + const int rows = m_qr.rows(); + ei_assert(b.rows() == rows); + result->resize(rows, b.cols()); + + // TODO(keir): There is almost certainly a faster way to multiply by + // 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) + .template marked() + .solveTriangularInPlace(result->corner(TopLeft, m_rank, result->cols())); + + return true; } /** \returns the matrix Q */ template -MatrixType QR::matrixQ(void) const +MatrixType QR::matrixQ() const { + ei_assert(m_isInitialized && "QR 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), ...] diff --git a/Eigen/src/QR/SelfAdjointEigenSolver.h b/Eigen/src/QR/SelfAdjointEigenSolver.h index 1c9ad513d..70984efab 100644 --- a/Eigen/src/QR/SelfAdjointEigenSolver.h +++ b/Eigen/src/QR/SelfAdjointEigenSolver.h @@ -1,5 +1,5 @@ // This file is part of Eigen, a lightweight C++ template library -// for linear algebra. Eigen itself is part of the KDE project. +// for linear algebra. // // Copyright (C) 2008 Gael Guennebaud // @@ -52,8 +52,8 @@ template class SelfAdjointEigenSolver typedef Tridiagonalization TridiagonalizationType; SelfAdjointEigenSolver() - : m_eivec(Size, Size), - m_eivalues(Size) + : m_eivec(int(Size), int(Size)), + m_eivalues(int(Size)) { ei_assert(Size!=Dynamic); } @@ -189,6 +189,14 @@ void SelfAdjointEigenSolver::compute(const MatrixType& matrix, bool assert(matrix.cols() == matrix.rows()); int n = matrix.cols(); m_eivalues.resize(n,1); + + if(n==1) + { + m_eivalues.coeffRef(0,0) = ei_real(matrix.coeff(0,0)); + m_eivec.setOnes(); + return; + } + m_eivec = matrix; // FIXME, should tridiag be a local variable of this function or an attribute of SelfAdjointEigenSolver ?