merge with tip

Eigen/CMakeLists.txt

@@ -1,4 +1,4 @@
-set(Eigen_HEADERS Core LU Cholesky QR Geometry Sparse Array SVD LeastSquares QtAlignedMalloc StdVector)
+set(Eigen_HEADERS Core LU Cholesky QR Geometry Sparse Array SVD LeastSquares QtAlignedMalloc StdVector Householder Jacobi)
 
 if(EIGEN_BUILD_LIB)
     set(Eigen_SRCS

Eigen/QR

@@ -35,7 +35,7 @@ namespace Eigen {
   * \endcode
   */
 
-#include "src/QR/QR.h"
+#include "src/QR/HouseholderQR.h"
 #include "src/QR/FullPivotingHouseholderQR.h"
 #include "src/QR/ColPivotingHouseholderQR.h"
 #include "src/QR/Tridiagonalization.h"

Eigen/src/Array/CwiseOperators.h

@@ -43,38 +43,6 @@ Cwise<ExpressionType>::sqrt() const
   return _expression();
 }
 
-/** \array_module
-  * 
-  * \returns an expression of the coefficient-wise exponential of *this.
-  *
-  * Example: \include Cwise_exp.cpp
-  * Output: \verbinclude Cwise_exp.out
-  *
-  * \sa pow(), log(), sin(), cos()
-  */
-template<typename ExpressionType>
-inline const EIGEN_CWISE_UNOP_RETURN_TYPE(ei_scalar_exp_op)
-Cwise<ExpressionType>::exp() const
-{
-  return _expression();
-}
-
-/** \array_module
-  * 
-  * \returns an expression of the coefficient-wise logarithm of *this.
-  *
-  * Example: \include Cwise_log.cpp
-  * Output: \verbinclude Cwise_log.out
-  *
-  * \sa exp()
-  */
-template<typename ExpressionType>
-inline const EIGEN_CWISE_UNOP_RETURN_TYPE(ei_scalar_log_op)
-Cwise<ExpressionType>::log() const
-{
-  return _expression();
-}
-
 /** \array_module
   * 
   * \returns an expression of the coefficient-wise cosine of *this.

Eigen/src/Array/Functors.h

@@ -69,40 +69,6 @@ struct ei_functor_traits<ei_scalar_sqrt_op<Scalar> >
   };
 };
 
-/** \internal
-  *
-  * \array_module
-  *
-  * \brief Template functor to compute the exponential of a scalar
-  *
-  * \sa class CwiseUnaryOp, Cwise::exp()
-  */
-template<typename Scalar> struct ei_scalar_exp_op EIGEN_EMPTY_STRUCT {
-  inline const Scalar operator() (const Scalar& a) const { return ei_exp(a); }
-  typedef typename ei_packet_traits<Scalar>::type Packet;
-  inline Packet packetOp(const Packet& a) const { return ei_pexp(a); }
-};
-template<typename Scalar>
-struct ei_functor_traits<ei_scalar_exp_op<Scalar> >
-{ enum { Cost = 5 * NumTraits<Scalar>::MulCost, PacketAccess = ei_packet_traits<Scalar>::HasExp }; };
-
-/** \internal
-  *
-  * \array_module
-  *
-  * \brief Template functor to compute the logarithm of a scalar
-  *
-  * \sa class CwiseUnaryOp, Cwise::log()
-  */
-template<typename Scalar> struct ei_scalar_log_op EIGEN_EMPTY_STRUCT {
-  inline const Scalar operator() (const Scalar& a) const { return ei_log(a); }
-  typedef typename ei_packet_traits<Scalar>::type Packet;
-  inline Packet packetOp(const Packet& a) const { return ei_plog(a); }
-};
-template<typename Scalar>
-struct ei_functor_traits<ei_scalar_log_op<Scalar> >
-{ enum { Cost = 5 * NumTraits<Scalar>::MulCost, PacketAccess = ei_packet_traits<Scalar>::HasLog }; };
-
 /** \internal
   *
   * \array_module

Eigen/src/CMakeLists.txt

@@ -7,3 +7,5 @@ ADD_SUBDIRECTORY(Array)
 ADD_SUBDIRECTORY(Geometry)
 ADD_SUBDIRECTORY(LeastSquares)
 ADD_SUBDIRECTORY(Sparse)
+ADD_SUBDIRECTORY(Jacobi)
+ADD_SUBDIRECTORY(Householder)

Eigen/src/Core/CwiseUnaryOp.h

@@ -205,6 +205,35 @@ MatrixBase<Derived>::cast() const
   return derived();
 }
 
+/** \returns an expression of the coefficient-wise exponential of *this.
+  *
+  * Example: \include Cwise_exp.cpp
+  * Output: \verbinclude Cwise_exp.out
+  *
+  * \sa pow(), log(), sin(), cos()
+  */
+template<typename ExpressionType>
+inline const EIGEN_CWISE_UNOP_RETURN_TYPE(ei_scalar_exp_op)
+Cwise<ExpressionType>::exp() const
+{
+  return _expression();
+}
+
+/** \returns an expression of the coefficient-wise logarithm of *this.
+  *
+  * Example: \include Cwise_log.cpp
+  * Output: \verbinclude Cwise_log.out
+  *
+  * \sa exp()
+  */
+template<typename ExpressionType>
+inline const EIGEN_CWISE_UNOP_RETURN_TYPE(ei_scalar_log_op)
+Cwise<ExpressionType>::log() const
+{
+  return _expression();
+}
+
+
 /** \relates MatrixBase */
 template<typename Derived>
 EIGEN_STRONG_INLINE const typename MatrixBase<Derived>::ScalarMultipleReturnType

Eigen/src/Core/Functors.h

@@ -298,6 +298,36 @@ template<typename Scalar>
 struct ei_functor_traits<ei_scalar_imag_op<Scalar> >
 { enum { Cost = 0, PacketAccess = false }; };
 
+/** \internal
+  *
+  * \brief Template functor to compute the exponential of a scalar
+  *
+  * \sa class CwiseUnaryOp, Cwise::exp()
+  */
+template<typename Scalar> struct ei_scalar_exp_op EIGEN_EMPTY_STRUCT {
+  inline const Scalar operator() (const Scalar& a) const { return ei_exp(a); }
+  typedef typename ei_packet_traits<Scalar>::type Packet;
+  inline Packet packetOp(const Packet& a) const { return ei_pexp(a); }
+};
+template<typename Scalar>
+struct ei_functor_traits<ei_scalar_exp_op<Scalar> >
+{ enum { Cost = 5 * NumTraits<Scalar>::MulCost, PacketAccess = ei_packet_traits<Scalar>::HasExp }; };
+
+/** \internal
+  *
+  * \brief Template functor to compute the logarithm of a scalar
+  *
+  * \sa class CwiseUnaryOp, Cwise::log()
+  */
+template<typename Scalar> struct ei_scalar_log_op EIGEN_EMPTY_STRUCT {
+  inline const Scalar operator() (const Scalar& a) const { return ei_log(a); }
+  typedef typename ei_packet_traits<Scalar>::type Packet;
+  inline Packet packetOp(const Packet& a) const { return ei_plog(a); }
+};
+template<typename Scalar>
+struct ei_functor_traits<ei_scalar_log_op<Scalar> >
+{ enum { Cost = 5 * NumTraits<Scalar>::MulCost, PacketAccess = ei_packet_traits<Scalar>::HasLog }; };
+
 /** \internal
   * \brief Template functor to multiply a scalar by a fixed other one
   *

Eigen/src/Jacobi/Jacobi.h

@@ -63,7 +63,7 @@ inline void MatrixBase<Derived>::applyJacobiOnTheRight(int p, int q, Scalar c, S
 /** Computes the cosine-sine pair (\a c, \a s) such that its associated
   * rotation \f$ J = ( \begin{array}{cc} c & s \\ -s' c \end{array} )\f$
   * applied to both the right and left of the 2x2 matrix
-  * \f$ B = ( \begin{array}{cc} x & y \\ _ & z \end{array} )\f$ yields
+  * \f$ B = ( \begin{array}{cc} x & y \\ * & z \end{array} )\f$ yields
   * a diagonal matrix A: \f$ A = J' B J \f$
   */
 template<typename Scalar>
@@ -149,7 +149,7 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY&
 
     const Packet pc = ei_pset1(c);
     const Packet ps = ei_pset1(s);
-    ei_conj_helper<true,false> cj;
+    ei_conj_helper<NumTraits<Scalar>::IsComplex,false> cj;
 
     for(int i=0; i<alignedStart; ++i)
     {

Eigen/src/LU/LU.h

@@ -534,7 +534,16 @@ bool LU<MatrixType>::solve(
 ) const
 {
   ei_assert(m_originalMatrix != 0 && "LU is not initialized.");
-  if(m_rank==0) return false;
+  result->resize(m_lu.cols(), b.cols());
+  if(m_rank==0)
+  {
+    if(b.squaredNorm() == RealScalar(0))
+    {
+      result->setZero();
+      return true;
+    }
+    else return false;
+  }
   
   /* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1}.
    * So we proceed as follows:
@@ -577,7 +586,6 @@ bool LU<MatrixType>::solve(
       .solveInPlace(c.corner(TopLeft, m_rank, c.cols()));
 
   // Step 4
-  result->resize(m_lu.cols(), b.cols());
   for(int i = 0; i < m_rank; ++i) result->row(m_q.coeff(i)) = c.row(i);
   for(int i = m_rank; i < m_lu.cols(); ++i) result->row(m_q.coeff(i)).setZero();
   return true;

Eigen/src/QR/ColPivotingHouseholderQR.h

@@ -31,14 +31,14 @@
   *
   * \class ColPivotingHouseholderQR
   *
-  * \brief Householder rank-revealing QR decomposition of a matrix
+  * \brief Householder rank-revealing QR decomposition of a matrix with column-pivoting
   *
   * \param MatrixType the type of the matrix of which we are computing the QR decomposition
   *
   * This class performs a rank-revealing QR decomposition using Householder transformations.
   *
-  * This decomposition performs full-pivoting in order to be rank-revealing and achieve optimal
+  * This decomposition performs column pivoting in order to be rank-revealing and improve
-  * numerical stability.
+  * numerical stability. It is slower than HouseholderQR, and faster than FullPivotingHouseholderQR.
   *
   * \sa MatrixBase::colPivotingHouseholderQr()
   */
@@ -82,6 +82,8 @@ template<typename MatrixType> class ColPivotingHouseholderQR
     /** 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.
       *
+      * \returns \c true if a solution exists, \c false if no solution 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.
@@ -95,13 +97,17 @@ template<typename MatrixType> class ColPivotingHouseholderQR
       * Output: \verbinclude ColPivotingHouseholderQR_solve.out
       */
     template<typename OtherDerived, typename ResultType>
-    void solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
+    bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
 
-    MatrixType matrixQ(void) const;
+    MatrixQType matrixQ(void) const;
 
     /** \returns a reference to the matrix where the Householder QR decomposition is stored
       */
-    const MatrixType& matrixQR() const { return m_qr; }
+    const MatrixType& matrixQR() const
+    {
+      ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized.");
+      return m_qr;
+    }
 
     ColPivotingHouseholderQR& compute(const MatrixType& matrix);
     
@@ -111,12 +117,122 @@ template<typename MatrixType> class ColPivotingHouseholderQR
       return m_cols_permutation;
     }
     
+    /** \returns the absolute value of the determinant of the matrix of which
+      * *this is the QR decomposition. It has only linear complexity
+      * (that is, O(n) where n is the dimension of the square matrix)
+      * as the QR decomposition has already been computed.
+      *
+      * \note This is only for square matrices.
+      *
+      * \warning a determinant can be very big or small, so for matrices
+      * of large enough dimension, there is a risk of overflow/underflow.
+      * One way to work around that is to use logAbsDeterminant() instead.
+      *
+      * \sa logAbsDeterminant(), MatrixBase::determinant()
+      */
+    typename MatrixType::RealScalar absDeterminant() const;
+
+    /** \returns the natural log of the absolute value of the determinant of the matrix of which
+      * *this is the QR decomposition. It has only linear complexity
+      * (that is, O(n) where n is the dimension of the square matrix)
+      * as the QR decomposition has already been computed.
+      *
+      * \note This is only for square matrices.
+      *
+      * \note This method is useful to work around the risk of overflow/underflow that's inherent
+      * to determinant computation.
+      *
+      * \sa absDeterminant(), MatrixBase::determinant()
+      */
+    typename MatrixType::RealScalar logAbsDeterminant() const;
+    
+    /** \returns the rank of the matrix of which *this is the QR decomposition.
+      *
+      * \note This is computed at the time of the construction of the QR decomposition. This
+      *       method does not perform any further computation.
+      */
     inline int rank() const
     {
       ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized.");
       return m_rank;
     }
 
+    /** \returns the dimension of the kernel of the matrix of which *this is the QR decomposition.
+      *
+      * \note Since the rank is computed at the time of the construction of the QR decomposition, this
+      *       method almost does not perform any further computation.
+      */
+    inline int dimensionOfKernel() const
+    {
+      ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized.");
+      return m_qr.cols() - m_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 at the time of the construction of the QR decomposition, this
+      *       method almost does not perform any further computation.
+      */
+    inline bool isInjective() const
+    {
+      ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized.");
+      return m_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 at the time of the construction of the QR decomposition, this
+      *       method almost does not perform any further computation.
+      */
+    inline bool isSurjective() const
+    {
+      ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized.");
+      return m_rank == m_qr.rows();
+    }
+
+    /** \returns true if the matrix of which *this is the QR decomposition is invertible.
+      *
+      * \note Since the rank is computed at the time of the construction of the QR decomposition, this
+      *       method almost does not perform any further computation.
+      */
+    inline bool isInvertible() const
+    {
+      ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized.");
+      return isInjective() && isSurjective();
+    }
+
+    /** Computes the inverse of the matrix of which *this is the QR decomposition.
+      *
+      * \param result a pointer to the matrix into which to store the inverse. Resized if needed.
+      *
+      * \note If this matrix is not invertible, *result is left with undefined coefficients.
+      *       Use isInvertible() to first determine whether this matrix is invertible.
+      *
+      * \sa inverse()
+      */
+    inline void computeInverse(MatrixType *result) const
+    {
+      ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized.");
+      ei_assert(m_qr.rows() == m_qr.cols() && "You can't take the inverse of a non-square matrix!");
+      solve(MatrixType::Identity(m_qr.rows(), m_qr.cols()), result);
+    }
+
+    /** \returns the inverse of the matrix of which *this is the QR decomposition.
+      *
+      * \note If this matrix is not invertible, the returned matrix has undefined coefficients.
+      *       Use isInvertible() to first determine whether this matrix is invertible.
+      *
+      * \sa computeInverse()
+      */
+    inline MatrixType inverse() const
+    {
+      MatrixType result;
+      computeInverse(&result);
+      return result;
+    }
+
   protected:
     MatrixType m_qr;
     HCoeffsType m_hCoeffs;
@@ -129,6 +245,22 @@ template<typename MatrixType> class ColPivotingHouseholderQR
 
 #ifndef EIGEN_HIDE_HEAVY_CODE
 
+template<typename MatrixType>
+typename MatrixType::RealScalar ColPivotingHouseholderQR<MatrixType>::absDeterminant() const
+{
+  ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized.");
+  ei_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
+  return ei_abs(m_qr.diagonal().prod());
+}
+
+template<typename MatrixType>
+typename MatrixType::RealScalar ColPivotingHouseholderQR<MatrixType>::logAbsDeterminant() const
+{
+  ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized.");
+  ei_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
+  return m_qr.diagonal().cwise().abs().cwise().log().sum();
+}
+
 template<typename MatrixType>
 ColPivotingHouseholderQR<MatrixType>& ColPivotingHouseholderQR<MatrixType>::compute(const MatrixType& matrix)
 {
@@ -199,12 +331,23 @@ ColPivotingHouseholderQR<MatrixType>& ColPivotingHouseholderQR<MatrixType>::comp
 
 template<typename MatrixType>
 template<typename OtherDerived, typename ResultType>
-void ColPivotingHouseholderQR<MatrixType>::solve(
+bool ColPivotingHouseholderQR<MatrixType>::solve(
   const MatrixBase<OtherDerived>& b,
   ResultType *result
 ) const
 {
   ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized.");
+  result->resize(m_qr.cols(), b.cols());
+  if(m_rank==0)
+  {
+    if(b.squaredNorm() == RealScalar(0))
+    {
+      result->setZero();
+      return true;
+    }
+    else return false;
+  }
+
   const int rows = m_qr.rows();
   const int cols = b.cols();
   ei_assert(b.rows() == rows);
@@ -219,18 +362,27 @@ void ColPivotingHouseholderQR<MatrixType>::solve(
      .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), m_hCoeffs.coeff(k), &temp.coeffRef(0));
   }
 
+  if(!isSurjective())
+  {
+    // is c is in the image of R ?
+    RealScalar biggest_in_upper_part_of_c = c.corner(TopLeft, m_rank, c.cols()).cwise().abs().maxCoeff();
+    RealScalar biggest_in_lower_part_of_c = c.corner(BottomLeft, rows-m_rank, c.cols()).cwise().abs().maxCoeff();
+    if(!ei_isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision*4))
+      return false;
+  }
+
   m_qr.corner(TopLeft, m_rank, m_rank)
       .template triangularView<UpperTriangular>()
       .solveInPlace(c.corner(TopLeft, m_rank, c.cols()));
 
-  result->resize(m_qr.cols(), b.cols());
   for(int i = 0; i < m_rank; ++i) result->row(m_cols_permutation.coeff(i)) = c.row(i);
   for(int i = m_rank; i < m_qr.cols(); ++i) result->row(m_cols_permutation.coeff(i)).setZero();
+  return true;
 }
 
 /** \returns the matrix Q */
 template<typename MatrixType>
-MatrixType ColPivotingHouseholderQR<MatrixType>::matrixQ() const
+typename ColPivotingHouseholderQR<MatrixType>::MatrixQType ColPivotingHouseholderQR<MatrixType>::matrixQ() const
 {
   ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized.");
   // compute the product H'_0 H'_1 ... H'_n-1,
@@ -239,7 +391,7 @@ MatrixType ColPivotingHouseholderQR<MatrixType>::matrixQ() const
   int rows = m_qr.rows();
   int cols = m_qr.cols();
   int size = std::min(rows,cols);
-  MatrixType res = MatrixType::Identity(rows, rows);
+  MatrixQType res = MatrixQType::Identity(rows, rows);
   Matrix<Scalar,1,MatrixType::RowsAtCompileTime> temp(rows);
   for (int k = size-1; k >= 0; k--)
   {

Eigen/src/QR/FullPivotingHouseholderQR.h

@@ -31,16 +31,16 @@
   *
   * \class FullPivotingHouseholderQR
   *
-  * \brief Householder rank-revealing QR decomposition of a matrix
+  * \brief Householder rank-revealing QR decomposition of a matrix with full pivoting
   *
   * \param MatrixType the type of the matrix of which we are computing the QR decomposition
   *
   * This class performs a rank-revealing QR decomposition using Householder transformations.
   *
-  * This decomposition performs full-pivoting in order to be rank-revealing and achieve optimal
+  * This decomposition performs a very prudent full pivoting in order to be rank-revealing and achieve optimal
-  * numerical stability.
+  * numerical stability. The trade-off is that it is slower than HouseholderQR and ColPivotingHouseholderQR.
   *
-  * \sa MatrixBase::householderRrqr()
+  * \sa MatrixBase::fullPivotingHouseholderQr()
   */
 template<typename MatrixType> class FullPivotingHouseholderQR
 {
@@ -62,12 +62,11 @@ template<typename MatrixType> class FullPivotingHouseholderQR
     typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType;
     typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType;
 
-    /**
+    /** \brief Default Constructor.
-    * \brief Default Constructor.
+      *
-    *
+      * The default constructor is useful in cases in which the user intends to
-    * The default constructor is useful in cases in which the user intends to
+      * perform decompositions via FullPivotingHouseholderQR::compute(const MatrixType&).
-    * perform decompositions via FullPivotingHouseholderQR::compute(const MatrixType&).
+      */
-    */
     FullPivotingHouseholderQR() : m_qr(), m_hCoeffs(), m_isInitialized(false) {}
 
     FullPivotingHouseholderQR(const MatrixType& matrix)
@@ -81,6 +80,8 @@ template<typename MatrixType> class FullPivotingHouseholderQR
     /** 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.
       *
+      * \returns \c true if a solution exists, \c false if no solution 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.
@@ -96,11 +97,15 @@ template<typename MatrixType> class FullPivotingHouseholderQR
     template<typename OtherDerived, typename ResultType>
     bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
 
-    MatrixType matrixQ(void) const;
+    MatrixQType matrixQ(void) const;
 
     /** \returns a reference to the matrix where the Householder QR decomposition is stored
       */
-    const MatrixType& matrixQR() const { return m_qr; }
+    const MatrixType& matrixQR() const
+    {
+      ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized.");
+      return m_qr;
+    }
 
     FullPivotingHouseholderQR& compute(const MatrixType& matrix);
     
@@ -125,11 +130,26 @@ template<typename MatrixType> class FullPivotingHouseholderQR
       *
       * \warning a determinant can be very big or small, so for matrices
       * of large enough dimension, there is a risk of overflow/underflow.
+      * One way to work around that is to use logAbsDeterminant() instead.
       *
-      * \sa MatrixBase::determinant()
+      * \sa logAbsDeterminant(), MatrixBase::determinant()
       */
     typename MatrixType::RealScalar absDeterminant() const;
 
+    /** \returns the natural log of the absolute value of the determinant of the matrix of which
+      * *this is the QR decomposition. It has only linear complexity
+      * (that is, O(n) where n is the dimension of the square matrix)
+      * as the QR decomposition has already been computed.
+      *
+      * \note This is only for square matrices.
+      *
+      * \note This method is useful to work around the risk of overflow/underflow that's inherent
+      * to determinant computation.
+      *
+      * \sa absDeterminant(), MatrixBase::determinant()
+      */
+    typename MatrixType::RealScalar logAbsDeterminant() const;
+    
     /** \returns the rank of the matrix of which *this is the QR decomposition.
       *
       * \note This is computed at the time of the construction of the QR decomposition. This
@@ -238,6 +258,14 @@ typename MatrixType::RealScalar FullPivotingHouseholderQR<MatrixType>::absDeterm
   return ei_abs(m_qr.diagonal().prod());
 }
 
+template<typename MatrixType>
+typename MatrixType::RealScalar FullPivotingHouseholderQR<MatrixType>::logAbsDeterminant() const
+{
+  ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized.");
+  ei_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
+  return m_qr.diagonal().cwise().abs().cwise().log().sum();
+}
+
 template<typename MatrixType>
 FullPivotingHouseholderQR<MatrixType>& FullPivotingHouseholderQR<MatrixType>::compute(const MatrixType& matrix)
 {
@@ -322,7 +350,16 @@ bool FullPivotingHouseholderQR<MatrixType>::solve(
 ) const
 {
   ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized.");
-  if(m_rank==0) return false;
+  result->resize(m_qr.cols(), b.cols());
+  if(m_rank==0)
+  {
+    if(b.squaredNorm() == RealScalar(0))
+    {
+      result->setZero();
+      return true;
+    }
+    else return false;
+  }
   
   const int rows = m_qr.rows();
   const int cols = b.cols();
@@ -351,7 +388,6 @@ bool FullPivotingHouseholderQR<MatrixType>::solve(
       .template triangularView<UpperTriangular>()
       .solveInPlace(c.corner(TopLeft, m_rank, c.cols()));
 
-  result->resize(m_qr.cols(), b.cols());
   for(int i = 0; i < m_rank; ++i) result->row(m_cols_permutation.coeff(i)) = c.row(i);
   for(int i = m_rank; i < m_qr.cols(); ++i) result->row(m_cols_permutation.coeff(i)).setZero();
   return true;
@@ -359,7 +395,7 @@ bool FullPivotingHouseholderQR<MatrixType>::solve(
 
 /** \returns the matrix Q */
 template<typename MatrixType>
-MatrixType FullPivotingHouseholderQR<MatrixType>::matrixQ() const
+typename FullPivotingHouseholderQR<MatrixType>::MatrixQType FullPivotingHouseholderQR<MatrixType>::matrixQ() const
 {
   ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized.");
   // compute the product H'_0 H'_1 ... H'_n-1,
@@ -368,7 +404,7 @@ MatrixType FullPivotingHouseholderQR<MatrixType>::matrixQ() const
   int rows = m_qr.rows();
   int cols = m_qr.cols();
   int size = std::min(rows,cols);
-  MatrixType res = MatrixType::Identity(rows, rows);
+  MatrixQType res = MatrixQType::Identity(rows, rows);
   Matrix<Scalar,1,MatrixType::RowsAtCompileTime> temp(rows);
   for (int k = size-1; k >= 0; k--)
   {

Eigen/src/QR/HouseholderQR.h

@@ -2,6 +2,7 @@
 // for linear algebra.
 //
 // Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
+// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
 //
 // Eigen is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -38,6 +39,10 @@
   * stored in a compact way compatible with LAPACK.
   *
   * Note that no pivoting is performed. This is \b not a rank-revealing decomposition.
+  * If you want that feature, use FullPivotingHouseholderQR or ColPivotingHouseholderQR instead.
+  *
+  * This Householder QR decomposition is faster, but less numerically stable and less feature-full than
+  * FullPivotingHouseholderQR or ColPivotingHouseholderQR.
   *
   * \sa MatrixBase::householderQr()
   */
@@ -46,15 +51,17 @@ template<typename MatrixType> class HouseholderQR
   public:
 
     enum {
-      MinSizeAtCompileTime = EIGEN_ENUM_MIN(MatrixType::ColsAtCompileTime,MatrixType::RowsAtCompileTime)
+      RowsAtCompileTime = MatrixType::RowsAtCompileTime,
+      ColsAtCompileTime = MatrixType::ColsAtCompileTime,
+      Options = MatrixType::Options,
+      DiagSizeAtCompileTime = EIGEN_ENUM_MIN(ColsAtCompileTime,RowsAtCompileTime)
     };
     
     typedef typename MatrixType::Scalar Scalar;
     typedef typename MatrixType::RealScalar RealScalar;
-    typedef Block<MatrixType, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> MatrixRBlockType;
+    typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixQType;
-    typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> MatrixTypeR;
+    typedef Matrix<Scalar, DiagSizeAtCompileTime, 1> HCoeffsType;
-    typedef Matrix<Scalar, MinSizeAtCompileTime, 1> HCoeffsType;
+    typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType;
-    typedef Matrix<Scalar, 1, MatrixType::ColsAtCompileTime> RowVectorType;
 
     /**
     * \brief Default Constructor.
@@ -72,15 +79,6 @@ template<typename MatrixType> class HouseholderQR
       compute(matrix);
     }
 
-    /** \returns a read-only expression of the matrix R of the actual the QR decomposition */
-    const TriangularView<NestByValue<MatrixRBlockType>, UpperTriangular>
-    matrixR(void) const
-    {
-      ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
-      int cols = m_qr.cols();
-      return MatrixRBlockType(m_qr, 0, 0, cols, cols).nestByValue().template triangularView<UpperTriangular>();
-    }
-
     /** 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.
       *
@@ -99,15 +97,48 @@ template<typename MatrixType> class HouseholderQR
     template<typename OtherDerived, typename ResultType>
     void solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
 
-    MatrixType matrixQ(void) const;
+    MatrixQType 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; }
+    const MatrixType& matrixQR() const
+    {
+        ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
+        return m_qr;
+    }
 
     HouseholderQR& compute(const MatrixType& matrix);
 
+    /** \returns the absolute value of the determinant of the matrix of which
+      * *this is the QR decomposition. It has only linear complexity
+      * (that is, O(n) where n is the dimension of the square matrix)
+      * as the QR decomposition has already been computed.
+      *
+      * \note This is only for square matrices.
+      *
+      * \warning a determinant can be very big or small, so for matrices
+      * of large enough dimension, there is a risk of overflow/underflow.
+      * One way to work around that is to use logAbsDeterminant() instead.
+      *
+      * \sa logAbsDeterminant(), MatrixBase::determinant()
+      */
+    typename MatrixType::RealScalar absDeterminant() const;
+
+    /** \returns the natural log of the absolute value of the determinant of the matrix of which
+      * *this is the QR decomposition. It has only linear complexity
+      * (that is, O(n) where n is the dimension of the square matrix)
+      * as the QR decomposition has already been computed.
+      *
+      * \note This is only for square matrices.
+      *
+      * \note This method is useful to work around the risk of overflow/underflow that's inherent
+      * to determinant computation.
+      *
+      * \sa absDeterminant(), MatrixBase::determinant()
+      */
+    typename MatrixType::RealScalar logAbsDeterminant() const;
+
   protected:
     MatrixType m_qr;
     HCoeffsType m_hCoeffs;
@@ -116,6 +147,22 @@ template<typename MatrixType> class HouseholderQR
 
 #ifndef EIGEN_HIDE_HEAVY_CODE
 
+template<typename MatrixType>
+typename MatrixType::RealScalar HouseholderQR<MatrixType>::absDeterminant() const
+{
+  ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
+  ei_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
+  return ei_abs(m_qr.diagonal().prod());
+}
+
+template<typename MatrixType>
+typename MatrixType::RealScalar HouseholderQR<MatrixType>::logAbsDeterminant() const
+{
+  ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
+  ei_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
+  return m_qr.diagonal().cwise().abs().cwise().log().sum();
+}
+
 template<typename MatrixType>
 HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType& matrix)
 {
@@ -177,7 +224,7 @@ void HouseholderQR<MatrixType>::solve(
 
 /** \returns the matrix Q */
 template<typename MatrixType>
-MatrixType HouseholderQR<MatrixType>::matrixQ() const
+typename HouseholderQR<MatrixType>::MatrixQType HouseholderQR<MatrixType>::matrixQ() const
 {
   ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
   // compute the product H'_0 H'_1 ... H'_n-1,
@@ -185,13 +232,13 @@ MatrixType HouseholderQR<MatrixType>::matrixQ() const
   // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...]
   int rows = m_qr.rows();
   int cols = m_qr.cols();
-  MatrixType res = MatrixType::Identity(rows, cols);
+  int size = std::min(rows,cols);
-  Matrix<Scalar,1,MatrixType::ColsAtCompileTime> temp(cols);
+  MatrixQType res = MatrixQType::Identity(rows, rows);
-  for (int k = cols-1; k >= 0; k--)
+  Matrix<Scalar,1,MatrixType::RowsAtCompileTime> temp(rows);
+  for (int k = size-1; k >= 0; k--)
   {
-    int remainingSize = rows-k;
+    res.block(k, k, rows-k, rows-k)
-    res.corner(BottomRight, remainingSize, cols-k)
+       .applyHouseholderOnTheLeft(m_qr.col(k).end(rows-k-1), ei_conj(m_hCoeffs.coeff(k)), &temp.coeffRef(k));
-       .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), ei_conj(m_hCoeffs.coeff(k)), &temp.coeffRef(k));
   }
   return res;
 }

test/qr.cpp

@@ -27,7 +27,6 @@
 
 template<typename MatrixType> void qr(const MatrixType& m)
 {
-  /* this test covers the following files: QR.h */
   int rows = m.rows();
   int cols = m.cols();
 
@@ -37,8 +36,11 @@ template<typename MatrixType> void qr(const MatrixType& m)
 
   MatrixType a = MatrixType::Random(rows,cols);
   HouseholderQR<MatrixType> qrOfA(a);
-  VERIFY_IS_APPROX(a, qrOfA.matrixQ() * qrOfA.matrixR().toDense());
+  MatrixType r = qrOfA.matrixQR();
-  VERIFY_IS_NOT_APPROX(a+MatrixType::Identity(rows, cols), qrOfA.matrixQ() * qrOfA.matrixR().toDense());
+  // FIXME need better way to construct trapezoid
+  for(int i = 0; i < rows; i++) for(int j = 0; j < cols; j++) if(i>j) r(i,j) = Scalar(0);
+
+  VERIFY_IS_APPROX(a, qrOfA.matrixQ() * r);
 
   SquareMatrixType b = a.adjoint() * a;
 
@@ -57,8 +59,9 @@ template<typename MatrixType> void qr(const MatrixType& m)
 
 template<typename MatrixType> void qr_invertible()
 {
-  /* this test covers the following files: QR.h */
   typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
+  typedef typename MatrixType::Scalar Scalar;
+
   int size = ei_random<int>(10,50);
 
   MatrixType m1(size, size), m2(size, size), m3(size, size);
@@ -75,6 +78,16 @@ template<typename MatrixType> void qr_invertible()
   m3 = MatrixType::Random(size,size);
   qr.solve(m3, &m2);
   VERIFY_IS_APPROX(m3, m1*m2);
+  
+  // now construct a matrix with prescribed determinant
+  m1.setZero();
+  for(int i = 0; i < size; i++) m1(i,i) = ei_random<Scalar>();
+  RealScalar absdet = ei_abs(m1.diagonal().prod());
+  m3 = qr.matrixQ(); // get a unitary
+  m1 = m3 * m1 * m3;
+  qr.compute(m1);
+  VERIFY_IS_APPROX(absdet, qr.absDeterminant());
+  VERIFY_IS_APPROX(ei_log(absdet), qr.logAbsDeterminant());
 }
 
 template<typename MatrixType> void qr_verify_assert()
@@ -82,9 +95,11 @@ template<typename MatrixType> void qr_verify_assert()
   MatrixType tmp;
 
   HouseholderQR<MatrixType> qr;
-  VERIFY_RAISES_ASSERT(qr.matrixR())
+  VERIFY_RAISES_ASSERT(qr.matrixQR())
   VERIFY_RAISES_ASSERT(qr.solve(tmp,&tmp))
   VERIFY_RAISES_ASSERT(qr.matrixQ())
+  VERIFY_RAISES_ASSERT(qr.absDeterminant())
+  VERIFY_RAISES_ASSERT(qr.logAbsDeterminant())
 }
 
 void test_qr()

test/qr_colpivoting.cpp

@@ -28,7 +28,6 @@
 
 template<typename MatrixType> void qr()
 {
-  /* this test covers the following files: QR.h */
   int rows = ei_random<int>(20,200), cols = ei_random<int>(20,200), cols2 = ei_random<int>(20,200);
   int rank = ei_random<int>(1, std::min(rows, cols)-1);
 
@@ -39,6 +38,10 @@ template<typename MatrixType> void qr()
   createRandomMatrixOfRank(rank,rows,cols,m1);
   ColPivotingHouseholderQR<MatrixType> qr(m1);
   VERIFY_IS_APPROX(rank, qr.rank());
+  VERIFY(cols - qr.rank() == qr.dimensionOfKernel());
+  VERIFY(!qr.isInjective());
+  VERIFY(!qr.isInvertible());
+  VERIFY(!qr.isSurjective());
   
   MatrixType r = qr.matrixQR();
   // FIXME need better way to construct trapezoid
@@ -54,14 +57,17 @@ template<typename MatrixType> void qr()
   MatrixType m2 = MatrixType::Random(cols,cols2);
   MatrixType m3 = m1*m2;
   m2 = MatrixType::Random(cols,cols2);
-  qr.solve(m3, &m2);
+  VERIFY(qr.solve(m3, &m2));
   VERIFY_IS_APPROX(m3, m1*m2);
+  m3 = MatrixType::Random(rows,cols2);
+  VERIFY(!qr.solve(m3, &m2));
 }
 
 template<typename MatrixType> void qr_invertible()
 {
-  /* this test covers the following files: RRQR.h */
   typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
+  typedef typename MatrixType::Scalar Scalar;
+
   int size = ei_random<int>(10,50);
 
   MatrixType m1(size, size), m2(size, size), m3(size, size);
@@ -78,6 +84,16 @@ template<typename MatrixType> void qr_invertible()
   m3 = MatrixType::Random(size,size);
   qr.solve(m3, &m2);
   VERIFY_IS_APPROX(m3, m1*m2);
+
+  // now construct a matrix with prescribed determinant
+  m1.setZero();
+  for(int i = 0; i < size; i++) m1(i,i) = ei_random<Scalar>();
+  RealScalar absdet = ei_abs(m1.diagonal().prod());
+  m3 = qr.matrixQ(); // get a unitary
+  m1 = m3 * m1 * m3;
+  qr.compute(m1);
+  VERIFY_IS_APPROX(absdet, qr.absDeterminant());
+  VERIFY_IS_APPROX(ei_log(absdet), qr.logAbsDeterminant());
 }
 
 template<typename MatrixType> void qr_verify_assert()
@@ -85,9 +101,17 @@ template<typename MatrixType> void qr_verify_assert()
   MatrixType tmp;
 
   ColPivotingHouseholderQR<MatrixType> qr;
-  VERIFY_RAISES_ASSERT(qr.matrixR())
+  VERIFY_RAISES_ASSERT(qr.matrixQR())
   VERIFY_RAISES_ASSERT(qr.solve(tmp,&tmp))
   VERIFY_RAISES_ASSERT(qr.matrixQ())
+  VERIFY_RAISES_ASSERT(qr.dimensionOfKernel())
+  VERIFY_RAISES_ASSERT(qr.isInjective())
+  VERIFY_RAISES_ASSERT(qr.isSurjective())
+  VERIFY_RAISES_ASSERT(qr.isInvertible())
+  VERIFY_RAISES_ASSERT(qr.computeInverse(&tmp))
+  VERIFY_RAISES_ASSERT(qr.inverse())
+  VERIFY_RAISES_ASSERT(qr.absDeterminant())
+  VERIFY_RAISES_ASSERT(qr.logAbsDeterminant())
 }
 
 void test_qr_colpivoting()

test/qr_fullpivoting.cpp

@@ -28,7 +28,6 @@
 
 template<typename MatrixType> void qr()
 {
-  /* this test covers the following files: QR.h */
   int rows = ei_random<int>(20,200), cols = ei_random<int>(20,200), cols2 = ei_random<int>(20,200);
   int rank = ei_random<int>(1, std::min(rows, cols)-1);
 
@@ -44,7 +43,6 @@ template<typename MatrixType> void qr()
   VERIFY(!qr.isInvertible());
   VERIFY(!qr.isSurjective());
 
-  
   MatrixType r = qr.matrixQR();
   // FIXME need better way to construct trapezoid
   for(int i = 0; i < rows; i++) for(int j = 0; j < cols; j++) if(i>j) r(i,j) = Scalar(0);
@@ -99,6 +97,7 @@ template<typename MatrixType> void qr_invertible()
   m1 = m3 * m1 * m3;
   qr.compute(m1);
   VERIFY_IS_APPROX(absdet, qr.absDeterminant());
+  VERIFY_IS_APPROX(ei_log(absdet), qr.logAbsDeterminant());
 }
 
 template<typename MatrixType> void qr_verify_assert()
@@ -106,9 +105,17 @@ template<typename MatrixType> void qr_verify_assert()
   MatrixType tmp;
 
   FullPivotingHouseholderQR<MatrixType> qr;
-  VERIFY_RAISES_ASSERT(qr.matrixR())
+  VERIFY_RAISES_ASSERT(qr.matrixQR())
   VERIFY_RAISES_ASSERT(qr.solve(tmp,&tmp))
   VERIFY_RAISES_ASSERT(qr.matrixQ())
+  VERIFY_RAISES_ASSERT(qr.dimensionOfKernel())
+  VERIFY_RAISES_ASSERT(qr.isInjective())
+  VERIFY_RAISES_ASSERT(qr.isSurjective())
+  VERIFY_RAISES_ASSERT(qr.isInvertible())
+  VERIFY_RAISES_ASSERT(qr.computeInverse(&tmp))
+  VERIFY_RAISES_ASSERT(qr.inverse())
+  VERIFY_RAISES_ASSERT(qr.absDeterminant())
+  VERIFY_RAISES_ASSERT(qr.logAbsDeterminant())
 }
 
 void test_qr_fullpivoting()