[eigen] qr().solve() implementation

[Keir Mierle <mierle@xxxxxxxxx>]

Here is a (incomplete) implementation of qr().solve(b, &x).

There is something wrong with the generated docs; it seems to take the
results from LU_solve.cpp rather than QR_solve.cpp but I don't have
time to figure this out now.

It is probably wasteful in that it explicitly computes matrixQ() and
matrixR(). However, it does work.

Keir

Index: doc/snippets/compile_snippet.cpp.in
===================================================================
--- doc/snippets/compile_snippet.cpp.in	(revision 917187)
+++ doc/snippets/compile_snippet.cpp.in	(working copy)
@@ -1,6 +1,7 @@
 #include <Eigen/Core>
 #include <Eigen/Array>
 #include <Eigen/LU>
+#include <Eigen/QR>
 #include <Eigen/Cholesky>
 #include <Eigen/Geometry>
 
Index: doc/snippets/QR_solve.cpp
===================================================================
--- doc/snippets/QR_solve.cpp	(revision 0)
+++ doc/snippets/QR_solve.cpp	(revision 0)
@@ -0,0 +1,14 @@
+typedef Matrix<float,3,3> Matrix3x3;
+Matrix3x3 m = Matrix3x3::Random();
+Matrix3f y = Matrix3f::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is the matrix y:" << endl << y << endl;
+Matrix3f x;
+if(m.qr().solve(y, &x))
+{
+  assert(y.isApprox(m*x));
+  cout << "Here is a solution x to the equation mx=y:" << endl << x << endl;
+}
+else
+  cout << "The equation mx=y does not have any solution." << endl;
+
Index: Eigen/src/QR/QR.h
===================================================================
--- Eigen/src/QR/QR.h	(revision 917187)
+++ Eigen/src/QR/QR.h	(working copy)
@@ -69,6 +69,32 @@
       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. Must be a vector
+      *          with the same number of rows as A.
+      * \param result a pointer to the vector 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:
@@ -160,6 +186,25 @@
   }
 }
 
+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);
+  ei_assert(b.cols() == 1);  // TODO(keir): Support more than one column in b.
+
+  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 = matrixR().solveTriangular(matrixQ().transpose()*b);
+  // TODO(keir): Triangular solver can't fail; is that what's wanted here?
+  return true;
+}
+
 /** \returns the matrix Q */
 template<typename MatrixType>
 MatrixType QR<MatrixType>::matrixQ(void) const