Re: [eigen] Blocked QR with column pivoting.

[Rasmus Larsen <rmlarsen@xxxxxxxxxx>]

Attached is a patch to the existing ColPivHouseholderQR that updates it to use the numerically stable norm downdate formula due to Drmac and Bujanovic (see http://www.netlib.org/lapack/lawnspdf/lawn176.pdf), which avoids the need to directly re-compute the norm of the pivot column every time.. This makes pivoting & rank determination stable for graded matrices, and the patch contains new and stricter unit tests. The change comes with a modest slowdown for small matrices (e.g. ~18% slower for 64x64). I think it is worth it, but perhaps you have ideas for how to improve my implementation? I tried vectorizing the norm downdate step, but it required multiple extra scratch vectors and didn't improve speed much.

Best,

   Rasmus

On Thu, Jan 21, 2016 at 2:27 AM, Gael Guennebaud <gael.guennebaud@xxxxxxxxx> wrote:

This is a very long time wanted feature. Sadly, nobody started to work on it yet. Retrospectively, this "soon" was really optimistic. A tile-based implementation would be even better (better multi-threading), and completing this class with complete orthogonalization for minimal norm solving would be a nice bonus too.

Any volunteer? or sponsorship?

gael 

On Thu, Jan 21, 2016 at 12:42 AM, Rasmus Larsen <rmlarsen@xxxxxxxxxx> wrote:

It would be great to have an implementation of the LAPACK algorithm in xGEQP3 in Eigen. In the table here:

http://eigen.tuxfamily.org/dox/group__TopicLinearAlgebraDecompositions.html

there is a "teaser" saying that ColPivHouseholderQR will soon be blocked. Is somebody actually working on this, and is so, what's the timeline? 

Rasmus

# HG changeset patch
# User Rasmus Munk Larsen <rmlarsen@xxxxxxxxxx>
# Date 1454022446 28800
#      Thu Jan 28 15:07:26 2016 -0800
# Node ID 9da6c621d055e425605362c92c5dd1f7768e9681
# Parent  71d2e8a2808ba022ccb88215ae548b8bf3ec065a
Change Eigen's ColPivHouseholderQR to use the numerically stable norm downdate formula from http://www.netlib.org/lapack/lawnspdf/lawn176.pdf, which has been used in LAPACK's xGEQPF and xGEQP3 since 2006. With the old formula, the code chooses the wrong pivots and fails to correctly determine rank on graded matrices.

This change also adds additional checks for non-increasing diagonal in R11 to existing unit tests, and adds a new unit test with the Kahan matrix, which consistently fails for the original code.

Benchmark timings on Intel(R) Xeon(R) CPU E5-1650 v3 @ 3.50GHz. Code compiled with AVX & FMA. I just ran on square matrices of 3 difference sizes.

Benchmark               Time(ns)     CPU(ns) Iterations
-------------------------------------------------------
Before:
BM_EigencolPivQR/64        53677       53627      12890
BM_EigencolPivQR/512    15265408    15250784         46
BM_EigencolPivQR/4k  15403556228 15388788368          2

After (non-vectorized version):
Benchmark               Time(ns)     CPU(ns) Iterations  Degradation
--------------------------------------------------------------------
BM_EigencolPivQR/64        63736       63669      10844         18.5%
BM_EigencolPivQR/512    16052546    16037381         43          5.1%
BM_EigencolPivQR/4k  15149263620 15132025316          2         -2.0%

Performance-wise there seems to be a ~18.5% degradation for small (64x64) matrices, probably due to the cost of more O(min(m,n)^2) sqrt operations that are not needed for the unstable formula.

diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h
--- a/Eigen/src/QR/ColPivHouseholderQR.h
+++ b/Eigen/src/QR/ColPivHouseholderQR.h
@@ -6,17 +6,17 @@
 //
 // This Source Code Form is subject to the terms of the Mozilla
 // Public License v. 2.0. If a copy of the MPL was not distributed
 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
 #ifndef EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
 #define EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
 
-namespace Eigen { 
+namespace Eigen {
 
 namespace internal {
 template<typename _MatrixType> struct traits<ColPivHouseholderQR<_MatrixType> >
  : traits<_MatrixType>
 {
   enum { Flags = 0 };
 };
 
@@ -26,21 +26,21 @@ template<typename _MatrixType> struct tr
   *
   * \class ColPivHouseholderQR
   *
   * \brief Householder rank-revealing QR decomposition of a matrix with column-pivoting
   *
   * \tparam _MatrixType the type of the matrix of which we are computing the QR decomposition
   *
   * This class performs a rank-revealing QR decomposition of a matrix \b A into matrices \b P, \b Q and \b R
-  * such that 
+  * such that
   * \f[
   *  \mathbf{A} \, \mathbf{P} = \mathbf{Q} \, \mathbf{R}
   * \f]
-  * by using Householder transformations. Here, \b P is a permutation matrix, \b Q a unitary matrix and \b R an 
+  * by using Householder transformations. Here, \b P is a permutation matrix, \b Q a unitary matrix and \b R an
   * upper triangular matrix.
   *
   * This decomposition performs column pivoting in order to be rank-revealing and improve
   * numerical stability. It is slower than HouseholderQR, and faster than FullPivHouseholderQR.
   *
   * \sa MatrixBase::colPivHouseholderQr()
   */
 template<typename _MatrixType> class ColPivHouseholderQR
@@ -62,75 +62,75 @@ template<typename _MatrixType> class Col
     typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, Options, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType;
     typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
     typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
     typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType;
     typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
     typedef typename internal::plain_row_type<MatrixType, RealScalar>::type RealRowVectorType;
     typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename HCoeffsType::ConjugateReturnType>::type> HouseholderSequenceType;
     typedef typename MatrixType::PlainObject PlainObject;
-    
+
   private:
-    
+
     typedef typename PermutationType::StorageIndex PermIndexType;
-    
+
   public:
 
     /**
     * \brief Default Constructor.
     *
     * The default constructor is useful in cases in which the user intends to
     * perform decompositions via ColPivHouseholderQR::compute(const MatrixType&).
     */
     ColPivHouseholderQR()
       : m_qr(),
         m_hCoeffs(),
         m_colsPermutation(),
         m_colsTranspositions(),
         m_temp(),
-        m_colSqNorms(),
+        m_colNorms(),
         m_isInitialized(false),
         m_usePrescribedThreshold(false) {}
 
     /** \brief Default Constructor with memory preallocation
       *
       * Like the default constructor but with preallocation of the internal data
       * according to the specified problem \a size.
       * \sa ColPivHouseholderQR()
       */
     ColPivHouseholderQR(Index rows, Index cols)
       : m_qr(rows, cols),
         m_hCoeffs((std::min)(rows,cols)),
         m_colsPermutation(PermIndexType(cols)),
         m_colsTranspositions(cols),
         m_temp(cols),
-        m_colSqNorms(cols),
+        m_colNorms(cols),
         m_isInitialized(false),
         m_usePrescribedThreshold(false) {}
 
     /** \brief Constructs a QR factorization from a given matrix
       *
       * This constructor computes the QR factorization of the matrix \a matrix by calling
       * the method compute(). It is a short cut for:
-      * 
+      *
       * \code
       * ColPivHouseholderQR<MatrixType> qr(matrix.rows(), matrix.cols());
       * qr.compute(matrix);
       * \endcode
-      * 
+      *
       * \sa compute()
       */
     template<typename InputType>
     explicit ColPivHouseholderQR(const EigenBase<InputType>& matrix)
       : m_qr(matrix.rows(), matrix.cols()),
         m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
         m_colsPermutation(PermIndexType(matrix.cols())),
         m_colsTranspositions(matrix.cols()),
         m_temp(matrix.cols()),
-        m_colSqNorms(matrix.cols()),
+        m_colNorms(matrix.cols()),
         m_isInitialized(false),
         m_usePrescribedThreshold(false)
     {
       compute(matrix.derived());
     }
 
     /** 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.
@@ -155,42 +155,42 @@ template<typename _MatrixType> class Col
     {
       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
       return Solve<ColPivHouseholderQR, Rhs>(*this, b.derived());
     }
 
     HouseholderSequenceType householderQ() const;
     HouseholderSequenceType matrixQ() const
     {
-      return householderQ(); 
+      return householderQ();
     }
 
     /** \returns a reference to the matrix where the Householder QR decomposition is stored
       */
     const MatrixType& matrixQR() const
     {
       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
       return m_qr;
     }
-    
-    /** \returns a reference to the matrix where the result Householder QR is stored 
-     * \warning The strict lower part of this matrix contains internal values. 
+
+    /** \returns a reference to the matrix where the result Householder QR is stored
+     * \warning The strict lower part of this matrix contains internal values.
      * Only the upper triangular part should be referenced. To get it, use
      * \code matrixR().template triangularView<Upper>() \endcode
-     * For rank-deficient matrices, use 
-     * \code 
-     * matrixR().topLeftCorner(rank(), rank()).template triangularView<Upper>() 
+     * For rank-deficient matrices, use
+     * \code
+     * matrixR().topLeftCorner(rank(), rank()).template triangularView<Upper>()
      * \endcode
      */
     const MatrixType& matrixR() const
     {
       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
       return m_qr;
     }
-    
+
     template<typename InputType>
     ColPivHouseholderQR& compute(const EigenBase<InputType>& matrix);
 
     /** \returns a const reference to the column permutation matrix */
     const PermutationType& colsPermutation() const
     {
       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
       return m_colsPermutation;
@@ -300,19 +300,19 @@ template<typename _MatrixType> class Col
     inline const Inverse<ColPivHouseholderQR> inverse() const
     {
       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
       return Inverse<ColPivHouseholderQR>(*this);
     }
 
     inline Index rows() const { return m_qr.rows(); }
     inline Index cols() const { return m_qr.cols(); }
-    
+
     /** \returns a const reference to the vector of Householder coefficients used to represent the factor \c Q.
-      * 
+      *
       * For advanced uses only.
       */
     const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
 
     /** Allows to prescribe a threshold to be used by certain methods, such as rank(),
       * who need to determine when pivots are to be considered nonzero. This is not used for the
       * QR decomposition itself.
       *
@@ -375,50 +375,50 @@ template<typename _MatrixType> class Col
       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
       return m_nonzero_pivots;
     }
 
     /** \returns the absolute value of the biggest pivot, i.e. the biggest
       *          diagonal coefficient of R.
       */
     RealScalar maxPivot() const { return m_maxpivot; }
-    
+
     /** \brief Reports whether the QR factorization was succesful.
       *
-      * \note This function always returns \c Success. It is provided for compatibility 
+      * \note This function always returns \c Success. It is provided for compatibility
       * with other factorization routines.
-      * \returns \c Success 
+      * \returns \c Success
       */
     ComputationInfo info() const
     {
       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
       return Success;
     }
-    
+
     #ifndef EIGEN_PARSED_BY_DOXYGEN
     template<typename RhsType, typename DstType>
     EIGEN_DEVICE_FUNC
     void _solve_impl(const RhsType &rhs, DstType &dst) const;
     #endif
 
   protected:
-    
+
     static void check_template_parameters()
     {
       EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
     }
-    
+
     void computeInPlace();
-    
+
     MatrixType m_qr;
     HCoeffsType m_hCoeffs;
     PermutationType m_colsPermutation;
     IntRowVectorType m_colsTranspositions;
     RowVectorType m_temp;
-    RealRowVectorType m_colSqNorms;
+    RealRowVectorType m_colNorms;
     bool m_isInitialized, m_usePrescribedThreshold;
     RealScalar m_prescribedThreshold, m_maxpivot;
     Index m_nonzero_pivots;
     Index m_det_pq;
 };
 
 template<typename MatrixType>
 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::absDeterminant() const
@@ -443,77 +443,73 @@ typename MatrixType::RealScalar ColPivHo
   *
   * \sa class ColPivHouseholderQR, ColPivHouseholderQR(const MatrixType&)
   */
 template<typename MatrixType>
 template<typename InputType>
 ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const EigenBase<InputType>& matrix)
 {
   check_template_parameters();
-  
+
   // the column permutation is stored as int indices, so just to be sure:
   eigen_assert(matrix.cols()<=NumTraits<int>::highest());
 
   m_qr = matrix;
-  
+
   computeInPlace();
-  
+
   return *this;
 }
 
 template<typename MatrixType>
 void ColPivHouseholderQR<MatrixType>::computeInPlace()
 {
   using std::abs;
+
   Index rows = m_qr.rows();
   Index cols = m_qr.cols();
   Index size = m_qr.diagonalSize();
-  
+
   m_hCoeffs.resize(size);
 
   m_temp.resize(cols);
 
   m_colsTranspositions.resize(m_qr.cols());
   Index number_of_transpositions = 0;
 
-  m_colSqNorms.resize(cols);
-  for(Index k = 0; k < cols; ++k)
-    m_colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm();
+  m_colNorms.resize(cols);
+  for (Index k = 0; k < cols; ++k)
+    m_colNorms.coeffRef(k) = m_qr.col(k).norm();
+  RealRowVectorType colNormsMostRecentDirect(m_colNorms);
 
-  RealScalar threshold_helper = m_colSqNorms.maxCoeff() * numext::abs2(NumTraits<Scalar>::epsilon()) / RealScalar(rows);
+  RealScalar threshold_helper =  numext::abs2(m_colNorms.maxCoeff() * NumTraits<Scalar>::epsilon()) / RealScalar(rows);
+  RealScalar norm_downdate_threshold = numext::sqrt(NumTraits<Scalar>::epsilon());
 
   m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
   m_maxpivot = RealScalar(0);
 
   for(Index k = 0; k < size; ++k)
   {
-    // first, we look up in our table m_colSqNorms which column has the biggest squared norm
+    // first, we look up in our table m_colNorms which column has the biggest norm
     Index biggest_col_index;
-    RealScalar biggest_col_sq_norm = m_colSqNorms.tail(cols-k).maxCoeff(&biggest_col_index);
+    RealScalar biggest_col_sq_norm = numext::abs2(m_colNorms.tail(cols-k).maxCoeff(&biggest_col_index));
     biggest_col_index += k;
 
-    // since our table m_colSqNorms accumulates imprecision at every step, we must now recompute
-    // the actual squared norm of the selected column.
-    // Note that not doing so does result in solve() sometimes returning inf/nan values
-    // when running the unit test with 1000 repetitions.
-    biggest_col_sq_norm = m_qr.col(biggest_col_index).tail(rows-k).squaredNorm();
-
-    // we store that back into our table: it can't hurt to correct our table.
-    m_colSqNorms.coeffRef(biggest_col_index) = biggest_col_sq_norm;
-
     // Track the number of meaningful pivots but do not stop the decomposition to make
     // sure that the initial matrix is properly reproduced. See bug 941.
     if(m_nonzero_pivots==size && biggest_col_sq_norm < threshold_helper * RealScalar(rows-k))
       m_nonzero_pivots = k;
 
     // apply the transposition to the columns
     m_colsTranspositions.coeffRef(k) = biggest_col_index;
     if(k != biggest_col_index) {
       m_qr.col(k).swap(m_qr.col(biggest_col_index));
-      std::swap(m_colSqNorms.coeffRef(k), m_colSqNorms.coeffRef(biggest_col_index));
+      std::swap(m_colNorms.coeffRef(k), m_colNorms.coeffRef(biggest_col_index));
+      std::swap(colNormsMostRecentDirect.coeffRef(k),
+                colNormsMostRecentDirect.coeffRef(biggest_col_index));
       ++number_of_transpositions;
     }
 
     // generate the householder vector, store it below the diagonal
     RealScalar beta;
     m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
 
     // apply the householder transformation to the diagonal coefficient
@@ -521,18 +517,39 @@ void ColPivHouseholderQR<MatrixType>::co
 
     // remember the maximum absolute value of diagonal coefficients
     if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta);
 
     // apply the householder transformation
     m_qr.bottomRightCorner(rows-k, cols-k-1)
         .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
 
-    // update our table of squared norms of the columns
-    m_colSqNorms.tail(cols-k-1) -= m_qr.row(k).tail(cols-k-1).cwiseAbs2();
+    // update our table of norms of the columns
+    for (Index j = k + 1; j < cols; ++j) {
+      // The following implements the stable norm downgrade step discussed in
+      // http://www.netlib.org/lapack/lawnspdf/lawn176.pdf
+      // and used in LAPACK routines xGEQPF and xGEQP3.
+      // See lines 278-297 in http://www.netlib.org/lapack/explore-html/dc/df4/sgeqpf_8f_source.html
+      if (m_colNorms.coeffRef(j) != 0) {
+        RealScalar temp = abs(m_qr.coeffRef(k, j)) / m_colNorms.coeffRef(j);
+        temp = (RealScalar(1) + temp) * (RealScalar(1) - temp);
+        temp = temp < 0 ? 0 : temp;
+        RealScalar temp2 =
+            temp * numext::abs2(m_colNorms.coeffRef(j) /
+                                colNormsMostRecentDirect.coeffRef(j));
+        if (temp2 <= norm_downdate_threshold) {
+          // The updated norm has become to inaccurate so re-compute the column
+          // norm directly.
+          m_colNorms.coeffRef(j) = m_qr.col(j).tail(rows - k - 1).norm();
+          colNormsMostRecentDirect.coeffRef(j) = m_colNorms.coeffRef(j);
+        } else {
+          m_colNorms.coeffRef(j) *= numext::sqrt(temp);
+        }
+      }
+    }
   }
 
   m_colsPermutation.setIdentity(PermIndexType(cols));
   for(PermIndexType k = 0; k < size/*m_nonzero_pivots*/; ++k)
     m_colsPermutation.applyTranspositionOnTheRight(k, PermIndexType(m_colsTranspositions.coeff(k)));
 
   m_det_pq = (number_of_transpositions%2) ? -1 : 1;
   m_isInitialized = true;
@@ -573,17 +590,17 @@ void ColPivHouseholderQR<_MatrixType>::_
 namespace internal {
 
 template<typename DstXprType, typename MatrixType, typename Scalar>
 struct Assignment<DstXprType, Inverse<ColPivHouseholderQR<MatrixType> >, internal::assign_op<Scalar>, Dense2Dense, Scalar>
 {
   typedef ColPivHouseholderQR<MatrixType> QrType;
   typedef Inverse<QrType> SrcXprType;
   static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar> &)
-  {    
+  {
     dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
   }
 };
 
 } // end namespace internal
 
 /** \returns the matrix Q as a sequence of householder transformations.
   * You can extract the meaningful part only by using:
diff --git a/test/qr_colpivoting.cpp b/test/qr_colpivoting.cpp
--- a/test/qr_colpivoting.cpp
+++ b/test/qr_colpivoting.cpp
@@ -14,16 +14,17 @@
 template<typename MatrixType> void qr()
 {
   typedef typename MatrixType::Index Index;
 
   Index rows = internal::random<Index>(2,EIGEN_TEST_MAX_SIZE), cols = internal::random<Index>(2,EIGEN_TEST_MAX_SIZE), cols2 = internal::random<Index>(2,EIGEN_TEST_MAX_SIZE);
   Index rank = internal::random<Index>(1, (std::min)(rows, cols)-1);
 
   typedef typename MatrixType::Scalar Scalar;
+  typedef typename MatrixType::RealScalar RealScalar;
   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> MatrixQType;
   MatrixType m1;
   createRandomPIMatrixOfRank(rank,rows,cols,m1);
   ColPivHouseholderQR<MatrixType> qr(m1);
   VERIFY_IS_EQUAL(rank, qr.rank());
   VERIFY_IS_EQUAL(cols - qr.rank(), qr.dimensionOfKernel());
   VERIFY(!qr.isInjective());
   VERIFY(!qr.isInvertible());
@@ -31,27 +32,46 @@ template<typename MatrixType> void qr()
 
   MatrixQType q = qr.householderQ();
   VERIFY_IS_UNITARY(q);
 
   MatrixType r = qr.matrixQR().template triangularView<Upper>();
   MatrixType c = q * r * qr.colsPermutation().inverse();
   VERIFY_IS_APPROX(m1, c);
 
+  // Verify that the absolute value of the diagonal elements in R are
+  // non-increasing until they reach the singularity threshold.
+  RealScalar threshold =
+      std::sqrt(RealScalar(rows)) * (std::abs)(r(0, 0)) * NumTraits<Scalar>::epsilon();
+  for (Index i = 0; i < (std::min)(rows, cols) - 1; ++i) {
+    RealScalar x = (std::abs)(r(i, i));
+    RealScalar y = (std::abs)(r(i + 1, i + 1));
+    if (x < threshold && y < threshold) continue;
+    if (test_isApproxOrLessThan(x, y)) {
+      for (Index j = 0; j < (std::min)(rows, cols); ++j) {
+        std::cout << "i = " << j << ", |r_ii| = " << (std::abs)(r(j, j)) << std::endl;
+      }
+      std::cout << "Failure at i=" << i << ", rank=" << rank
+                << ", threshold=" << threshold << std::endl;
+    }
+    VERIFY_IS_APPROX_OR_LESS_THAN(y, x);
+  }
+
   MatrixType m2 = MatrixType::Random(cols,cols2);
   MatrixType m3 = m1*m2;
   m2 = MatrixType::Random(cols,cols2);
   m2 = qr.solve(m3);
   VERIFY_IS_APPROX(m3, m1*m2);
 }
 
 template<typename MatrixType, int Cols2> void qr_fixedsize()
 {
   enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime };
   typedef typename MatrixType::Scalar Scalar;
+  typedef typename MatrixType::RealScalar RealScalar;
   int rank = internal::random<int>(1, (std::min)(int(Rows), int(Cols))-1);
   Matrix<Scalar,Rows,Cols> m1;
   createRandomPIMatrixOfRank(rank,Rows,Cols,m1);
   ColPivHouseholderQR<Matrix<Scalar,Rows,Cols> > qr(m1);
   VERIFY_IS_EQUAL(rank, qr.rank());
   VERIFY_IS_EQUAL(Cols - qr.rank(), qr.dimensionOfKernel());
   VERIFY_IS_EQUAL(qr.isInjective(), (rank == Rows));
   VERIFY_IS_EQUAL(qr.isSurjective(), (rank == Cols));
@@ -61,16 +81,77 @@ template<typename MatrixType, int Cols2>
   Matrix<Scalar,Rows,Cols> c = qr.householderQ() * r * qr.colsPermutation().inverse();
   VERIFY_IS_APPROX(m1, c);
 
   Matrix<Scalar,Cols,Cols2> m2 = Matrix<Scalar,Cols,Cols2>::Random(Cols,Cols2);
   Matrix<Scalar,Rows,Cols2> m3 = m1*m2;
   m2 = Matrix<Scalar,Cols,Cols2>::Random(Cols,Cols2);
   m2 = qr.solve(m3);
   VERIFY_IS_APPROX(m3, m1*m2);
+  // Verify that the absolute value of the diagonal elements in R are
+  // non-increasing until they reache the singularity threshold.
+  RealScalar threshold =
+      std::sqrt(RealScalar(Rows)) * (std::abs)(r(0, 0)) * NumTraits<Scalar>::epsilon();
+  for (Index i = 0; i < (std::min)(int(Rows), int(Cols)) - 1; ++i) {
+    RealScalar x = (std::abs)(r(i, i));
+    RealScalar y = (std::abs)(r(i + 1, i + 1));
+    if (x < threshold && y < threshold) continue;
+    if (test_isApproxOrLessThan(x, y)) {
+      for (Index j = 0; j < (std::min)(int(Rows), int(Cols)); ++j) {
+        std::cout << "i = " << j << ", |r_ii| = " << (std::abs)(r(j, j)) << std::endl;
+      }
+      std::cout << "Failure at i=" << i << ", rank=" << rank
+                << ", threshold=" << threshold << std::endl;
+    }
+    VERIFY_IS_APPROX_OR_LESS_THAN(y, x);
+  }
+}
+
+// This test is meant to verify that pivots are chosen such that
+// even for a graded matrix, the diagonal of R falls of roughly
+// monotonically until it reaches the threshold for singularity.
+// We use the so-called Kahan matrix, which is a famous counter-example
+// for rank-revealing QR. See
+// http://www.netlib.org/lapack/lawnspdf/lawn176.pdf
+// page 3 for more detail.
+template<typename MatrixType> void qr_kahan_matrix()
+{
+  typedef typename MatrixType::Index Index;
+  typedef typename MatrixType::Scalar Scalar;
+  typedef typename MatrixType::RealScalar RealScalar;
+
+  Index rows = 300, cols = rows;
+
+  MatrixType m1;
+  m1.setZero(rows,cols);
+  RealScalar s = std::pow(NumTraits<RealScalar>::epsilon(), 1.0 / rows);
+  RealScalar c = std::sqrt(1 - s*s);
+  for (Index i = 0; i < rows; ++i) {
+    m1(i, i) = pow(s, i);
+    m1.row(i).tail(rows - i - 1) = -pow(s, i) * c * MatrixType::Ones(1, rows - i - 1);
+  }
+  m1 = (m1 + m1.transpose()).eval();
+  ColPivHouseholderQR<MatrixType> qr(m1);
+  MatrixType r = qr.matrixQR().template triangularView<Upper>();
+
+  RealScalar threshold =
+      std::sqrt(RealScalar(rows)) * (std::abs)(r(0, 0)) * NumTraits<Scalar>::epsilon();
+  for (Index i = 0; i < (std::min)(rows, cols) - 1; ++i) {
+    RealScalar x = (std::abs)(r(i, i));
+    RealScalar y = (std::abs)(r(i + 1, i + 1));
+    if (x < threshold && y < threshold) continue;
+    if (test_isApproxOrLessThan(x, y)) {
+      for (Index j = 0; j < (std::min)(rows, cols); ++j) {
+        std::cout << "i = " << j << ", |r_ii| = " << (std::abs)(r(j, j)) << std::endl;
+      }
+      std::cout << "Failure at i=" << i << ", rank=" << qr.rank()
+                << ", threshold=" << threshold << std::endl;
+    }
+    VERIFY_IS_APPROX_OR_LESS_THAN(y, x);
+  }
 }
 
 template<typename MatrixType> void qr_invertible()
 {
   using std::log;
   using std::abs;
   typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
   typedef typename MatrixType::Scalar Scalar;
@@ -127,16 +208,21 @@ void test_qr_colpivoting()
     CALL_SUBTEST_2( qr<MatrixXd>() );
     CALL_SUBTEST_3( qr<MatrixXcd>() );
     CALL_SUBTEST_4(( qr_fixedsize<Matrix<float,3,5>, 4 >() ));
     CALL_SUBTEST_5(( qr_fixedsize<Matrix<double,6,2>, 3 >() ));
     CALL_SUBTEST_5(( qr_fixedsize<Matrix<double,1,1>, 1 >() ));
   }
 
   for(int i = 0; i < g_repeat; i++) {
+    CALL_SUBTEST_1( qr_kahan_matrix<MatrixXf>() );
+    CALL_SUBTEST_2( qr_kahan_matrix<MatrixXd>() );
+  }
+
+  for(int i = 0; i < g_repeat; i++) {
     CALL_SUBTEST_1( qr_invertible<MatrixXf>() );
     CALL_SUBTEST_2( qr_invertible<MatrixXd>() );
     CALL_SUBTEST_6( qr_invertible<MatrixXcf>() );
     CALL_SUBTEST_3( qr_invertible<MatrixXcd>() );
   }
 
   CALL_SUBTEST_7(qr_verify_assert<Matrix3f>());
   CALL_SUBTEST_8(qr_verify_assert<Matrix3d>());