modernize HouseholderQR too, uniformize all that stuff, update tests

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/QR/ColPivotingHouseholderQR.h

@@ -99,11 +99,15 @@ template<typename MatrixType> class ColPivotingHouseholderQR
     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 && "ColPivotingHouseholderQR is not initialized.");
+      return m_qr;
+    }
 
     ColPivotingHouseholderQR& compute(const MatrixType& matrix);
     
@@ -363,9 +367,10 @@ bool ColPivotingHouseholderQR<MatrixType>::solve(
     // 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))
+    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()));
@@ -377,7 +382,7 @@ bool ColPivotingHouseholderQR<MatrixType>::solve(
 
 /** \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,
@@ -386,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

@@ -97,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);
     
@@ -391,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,
@@ -400,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

@@ -101,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

@@ -105,7 +105,7 @@ 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())
@@ -114,6 +114,8 @@ template<typename MatrixType> void qr_verify_assert()
   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()