[eigen] Re: Bug(s) report: SparseQR on tall-thin matrices

[Julian Kent <julian.kent@xxxxxxxxx>]

Just a follow up, I think I've fixed bug 2, the following patch should work:

+++ SparseQR.h  2017-01-10 09:31:56.600356224 +0100
@@ -609,7 +609,7 @@
   // Get the references
   SparseQR_QProduct(const SparseQRType& qr, const Derived& other, bool transpose) :
   m_qr(qr),m_other(other),m_transpose(transpose) {}
-  inline Index rows() const { return m_transpose ? m_qr.rows() : m_qr.cols(); }
+  inline Index rows() const { return m_transpose ? m_qr.cols() : m_qr.rows(); }
   inline Index cols() const { return m_other.cols(); }
  
   // Assign to a vector
@@ -619,7 +619,12 @@
     Index m = m_qr.rows();
     Index n = m_qr.cols();
     Index diagSize = (std::min)(m,n);
-    res = m_other;
+    res.topLeftCorner(m_other.rows(), m_other.cols()) = m_other;
+
+    //set the rest to 0
+    res.topRightCorner(m_other.rows(), res.cols() - m_other.cols()).fill(0);
+    res.bottomLeftCorner(res.rows() - m_other.rows(), m_other.cols()).fill(0);
+    res.bottomRightCorner(res.rows() - m_other.rows(), res.cols() - m_other.cols()).fill(0);
     if (m_transpose)
     {
       eigen_assert(m_qr.m_Q.rows() == m_other.rows() && "Non conforming object sizes");
@@ -637,7 +642,7 @@
     }
     else
     {
-      eigen_assert(m_qr.m_Q.rows() == m_other.rows() && "Non conforming object sizes");
+      eigen_assert(m_qr.m_Q.cols() == m_other.rows() && "Non conforming object sizes");
       // Compute res = Q * other column by column
       for(Index j = 0; j < res.cols(); j++)
       {

This solution should work when &res == &m_other as well, since the topRight/bottomLeft/bottomRight corner fills should go to no-ops, since one or both of the width/height of the block will be 0.
I'm not sure whether the repeated fill(0) calls are the best way to do this, basically we need to set everything not in the top left corner to 0. If there was some kind of invertSelection() it would be much easier to read.

For bug 1) I am just taking the matrixR().topRows(n) manually in other code. I'm not sure where would be best to put a conservativeResize(n, n) - at the end of factorize(...) maybe?

Anyway, combining these two fixes gives me the correct results.

Cheers

Julian

On Mon, Jan 9, 2017 at 3:51 PM, Julian Kent <julian.kent@xxxxxxxxx> wrote:

While trying to use SparseQR on a matrix A with rows > cols, I found 2 bugs:

1) The size of qr.matrixR() is m x n, instead of n x n as expected. SparseQR.h:305 initialises m_R with size (m,n), and nothing does any resizing. For now I'm just taking the topRows(n), but I'm not entirely sure this is correct, and it certainly isn't the behaviour I expect. Shouldn't there be a non-destructive resize, if the extra rows are really necessary for intermediate procesing?

2) qr.matrixQ() claims to be size m x n, as expected. However, trying to multiply qr.matrixQ() with a n x k dense matrix gives an assertion error:
Eigen/src/SparseQR/SparseQR.h:640: void Eigen::SparseQR_QProduct<SparseQRType, Derived>::evalTo(DesType&) const [with DesType = Eigen::Matrix<double, -1, -1>; SparseQRType = Eigen::SparseQR<Eigen::SparseMatrix<double>, Eigen::NaturalOrdering<int> >; Derived = Eigen::Matrix<double, -1, -1>]: Assertion `m_qr.m_Q.rows() == m_other.rows() && "Non conforming object sizes"' failed.

In the .solve(...) only matrixQ.transpose() is used, which is probably why this hasn't shown up earlier.

These bugs may be interacting with each other to fool any accuracy tests using A*P = Q*R on tall-thin matrices, with the extra rows in R passing the assert in Q.

Let me know if you need example matrices to work with.

Thanks

Julian Kent