QR: add isInjective(), isSurjective(),

    mark isFullRank() deprecated,
    add solve() (mix of Keir's patch and LU::solve())
=> there is big problem with complex which are not working

4 files changed, 225 insertions, 23 deletions

diff --git a/Eigen/src/QR/QR.h b/Eigen/src/QR/QR.h
index 13d49a4a3..0ea839d25 100644
--- a/Eigen/src/QR/QR.h
+++ b/Eigen/src/QR/QR.h
@@ -55,12 +55,64 @@ template<typename MatrixType> class QR
     {
       _compute(matrix);
     }
-
-    /** \returns whether or not the matrix is of full rank */
-    bool isFullRank() const { return rank() == std::min(m_qr.rows(),m_qr.cols()); }
     
+    /** \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.
+      */
+    bool isFullRank() const EIGEN_DEPRECATED { 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
+    {
+      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
+    {
+      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
+    {
+      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
+    {
+      return isInjective() && isSurjective();
+    }
+    
     /** \returns a read-only expression of the matrix R of the actual the QR decomposition */
     const Part<NestByValue<MatrixRBlockType>, UpperTriangular>
     matrixR(void) const
@@ -69,6 +121,32 @@ template<typename MatrixType> class QR
       return MatrixRBlockType(m_qr, 0, 0, cols, cols).nestByValue().template part<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.
+      *
+      * \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<typename OtherDerived, typename ResultType>
+    bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
+
     MatrixType matrixQ(void) const;
 
   private:
@@ -88,12 +166,11 @@ int QR<MatrixType>::rank() const
 {
   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<n; ++i)
-      if (ei_isMuchSmallerThan(m_qr.diagonal().coeff(i), maxCoeff))
-        --m_rank;
+    RealScalar maxCoeff = m_qr.diagonal().cwise().abs().maxCoeff();
+    int n = m_qr.cols();
+    m_rank = 0;
+    while(m_rank<n && !ei_isMuchSmallerThan(m_qr.diagonal().coeff(m_rank), maxCoeff))
+      ++m_rank;
     m_rankIsUptodate = true;
   }
   return m_rank;
@@ -132,7 +209,7 @@ void QR<MatrixType>::_compute(const MatrixType& matrix)
         m_hCoeffs.coeffRef(k) = 0;
       }
     }
-    else if ( (!ei_isMuchSmallerThan(beta=m_qr.col(k).end(remainingSize-1).squaredNorm(),static_cast<Scalar>(1))) || ei_imag(v0)==0 )
+    else if ( (!ei_isMuchSmallerThan(beta=m_qr.col(k).end(remainingSize-1).squaredNorm(),static_cast<Scalar>(1))) ) // FIXME what about ei_imag(v0) ??
     {
       // form k-th Householder vector
       beta = ei_sqrt(ei_abs2(v0)+beta);
@@ -160,9 +237,40 @@ void QR<MatrixType>::_compute(const MatrixType& matrix)
   }
 }
 
+template<typename MatrixType>
+template<typename OtherDerived, typename ResultType>
+bool QR<MatrixType>::solve(
+  const MatrixBase<OtherDerived>& b,
+  ResultType *result
+) const
+{
+  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<UpperTriangular>()
+      .solveTriangularInPlace(result->corner(TopLeft, m_rank, result->cols()));
+
+  return true;
+}
+
 /** \returns the matrix Q */
 template<typename MatrixType>
-MatrixType QR<MatrixType>::matrixQ(void) const
+MatrixType QR<MatrixType>::matrixQ() const
 {
   // 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'
diff --git a/doc/snippets/compile_snippet.cpp.in b/doc/snippets/compile_snippet.cpp.in
index 3eaee98ac..d074cac50 100644
--- a/doc/snippets/compile_snippet.cpp.in
+++ b/doc/snippets/compile_snippet.cpp.in
@@ -1,6 +1,7 @@
 #include <Eigen/Core>
 #include <Eigen/Array>
 #include <Eigen/LU>
+#include <Eigen/QR>
 #include <Eigen/Cholesky>
 #include <Eigen/Geometry>
 
diff --git a/test/qr.cpp b/test/qr.cpp
index 9ce808429..2b9579c9c 100644
--- a/test/qr.cpp
+++ b/test/qr.cpp
@@ -27,9 +27,7 @@
 
 template<typename MatrixType> void qr(const MatrixType& m)
 {
-  /* this test covers the following files:
-     QR.h
-  */
+  /* this test covers the following files: QR.h */
   int rows = m.rows();
   int cols = m.cols();
 
@@ -57,6 +55,98 @@ template<typename MatrixType> void qr(const MatrixType& m)
   VERIFY_IS_APPROX(b, hess.matrixQ() * hess.matrixH() * hess.matrixQ().adjoint());
 }
 
+template<typename Derived>
+void doSomeRankPreservingOperations(Eigen::MatrixBase<Derived>& m)
+{
+  typedef typename Derived::RealScalar RealScalar;
+  for(int a = 0; a < 3*(m.rows()+m.cols()); a++)
+  {
+    RealScalar d = Eigen::ei_random<RealScalar>(-1,1);
+    int i = Eigen::ei_random<int>(0,m.rows()-1); // i is a random row number
+    int j;
+    do {
+      j = Eigen::ei_random<int>(0,m.rows()-1);
+    } while (i==j); // j is another one (must be different)
+    m.row(i) += d * m.row(j);
+
+    i = Eigen::ei_random<int>(0,m.cols()-1); // i is a random column number
+    do {
+      j = Eigen::ei_random<int>(0,m.cols()-1);
+    } while (i==j); // j is another one (must be different)
+    m.col(i) += d * m.col(j);
+  }
+}
+
+template<typename MatrixType> void qr_non_invertible()
+{
+  /* this test covers the following files: QR.h */
+  // NOTE there seems to be a problem with too small sizes -- could easily lie in the doSomeRankPreservingOperations function
+  int rows = ei_random<int>(20,200), cols = ei_random<int>(20,rows), cols2 = ei_random<int>(20,rows);
+  int rank = ei_random<int>(1, std::min(rows, cols)-1);
+
+  MatrixType m1(rows, cols), m2(cols, cols2), m3(rows, cols2), k(1,1);
+  m1 = MatrixType::Random(rows,cols);
+  if(rows <= cols)
+    for(int i = rank; i < rows; i++) m1.row(i).setZero();
+  else
+    for(int i = rank; i < cols; i++) m1.col(i).setZero();
+  doSomeRankPreservingOperations(m1);
+
+  QR<MatrixType> lu(m1);
+//   typename LU<MatrixType>::KernelResultType m1kernel = lu.kernel();
+//   typename LU<MatrixType>::ImageResultType m1image = lu.image();
+
+  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<typename MatrixType> void qr_invertible()
+{
+  /* this test covers the following files: QR.h */
+  typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
+  int size = ei_random<int>(10,200);
+
+  MatrixType m1(size, size), m2(size, size), m3(size, size);
+  m1 = MatrixType::Random(size,size);
+
+  if (ei_is_same_type<RealScalar,float>::ret)
+  {
+    // let's build a matrix more stable to inverse
+    MatrixType a = MatrixType::Random(size,size*2);
+    m1 += a * a.adjoint();
+  }
+
+  QR<MatrixType> lu(m1);
+  VERIFY(0 == lu.dimensionOfKernel());
+  VERIFY(size == lu.rank());
+  VERIFY(lu.isInjective());
+  VERIFY(lu.isSurjective());
+  VERIFY(lu.isInvertible());
+//   VERIFY(lu.image().lu().isInvertible());
+  m3 = MatrixType::Random(size,size);
+  lu.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));
+}
+
 void test_qr()
 {
   for(int i = 0; i < 1; i++) {
@@ -66,13 +156,17 @@ void test_qr()
     CALL_SUBTEST( qr(MatrixXcd(5,5)) );
     CALL_SUBTEST( qr(MatrixXcd(7,3)) );
   }
-
-  // small isFullRank test
-  {
-    Matrix3d mat;
-    mat << 1, 45, 1, 2, 2, 2, 1, 2, 3;
-    VERIFY(mat.qr().isFullRank());
-    mat << 1, 1, 1, 2, 2, 2, 1, 2, 3;
-    VERIFY(!mat.qr().isFullRank());
+  
+  for(int i = 0; i < g_repeat; i++) {
+    CALL_SUBTEST( qr_non_invertible<MatrixXf>() );
+    CALL_SUBTEST( qr_non_invertible<MatrixXd>() );
+    // TODO fix issue with complex
+//     CALL_SUBTEST( qr_non_invertible<MatrixXcf>() );
+//     CALL_SUBTEST( qr_non_invertible<MatrixXcd>() );
+    CALL_SUBTEST( qr_invertible<MatrixXf>() );
+    CALL_SUBTEST( qr_invertible<MatrixXd>() );
+    // TODO fix issue with complex
+//     CALL_SUBTEST( qr_invertible<MatrixXcf>() );
+//     CALL_SUBTEST( qr_invertible<MatrixXcd>() );
   }
 }
diff --git a/test/sparse_basic.cpp b/test/sparse_basic.cpp
index 08f19d7e4..23b8526b7 100644
--- a/test/sparse_basic.cpp
+++ b/test/sparse_basic.cpp
@@ -193,7 +193,6 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re
         }
       }
       m2.endFill();
-      //std::cerr << m1 << "\n\n" << m2 << "\n";
       VERIFY_IS_APPROX(m2,m1);
     }