New user self-registration is currently disabled. Please email eigen-core-team @ lists.tuxfamily.org if you need an account.
Bug 836 - Wrong results of SparseQR for non-square matrices
Wrong results of SparseQR for non-square matrices
Status: RESOLVED FIXED
Product: Eigen
Classification: Unclassified
Component: Sparse
3.2
All All
: Normal Wrong Result
Assigned To: Nobody
:
Depends on:
Blocks:
  Show dependency treegraph
 
Reported: 2014-06-28 03:02 UTC by Henrik Friberg
Modified: 2014-07-01 22:53 UTC (History)
2 users (show)



Attachments

Description Henrik Friberg 2014-06-28 03:02:19 UTC
I suspect the SparseQR module was never tested for non-square matrices.

Example 1:
Take the code at end of this bugreport with
A << 1, 0, 0, 0,
     0, 1, 0, 0,
     0, 0, 1, 1;

Not only is AP not equal to QR in code below, but I also get
QR.matrixQ().rows() = 3, 
QR.matrixQ().cols() = 4,
while,
Q.rows() = 3, 
Q.cols() = 3,
and
QR.matrixR().rows() = 4, 
QR.matrixR().cols() = 4.

This is the example should be easy to debug because there should essentially be nothing going on. The MATLAB solution is P = I [4x4], Q = I [3x3], R = A [3 x 4].

Example 2:
Take the code at end of this bugreport with
A << 1, 1, 0, 0,
     0, 1, 0, 0,
     0, 0, 1, 1;

The same thing happens as above, only this time the QR.rank() is 2 which is wrong unless it is another definition than I am use to.

------------------
SparseQR<SparseMatrix<double, ColMajor>, COLAMDOrdering<int> > QR;
MatrixXd A(3, 4);

A << 1, ?, 0, 0,
     0, 1, 0, 0,
     0, 0, 1, 1;

QR.compute(A.sparseView());

//
// QR rank
//
printf(", %d", QR.rank());

//
// AP matrix
//
printf("=[%d x %d]", AP.rows(), AP.cols());
for (i = 0; i < AP.outerSize(); ++i) {
  for (SparseMatrix<double>::InnerIterator it(AP, i); it; ++it) {
    printf("(%d %d %g)", it.row(), it.col(), it.value());
  }
}

// 
// Q matrix
//
printf("[%d x %d]", QR.matrixQ().rows(), QR.matrixQ().cols());
Q = QR.matrixQ();
printf("=[%d x %d]", Q.rows(), Q.cols());
for (i = 0; i < Q.outerSize(); ++i) {
  for (SparseMatrix<double>::InnerIterator it(Q, i); it; ++it) {
    printf("(%d %d %g)", it.row(), it.col(), it.value());
  }
}

//
// R matrix
//
printf("[%d x %d]", QR.matrixR().rows(), QR.matrixR().cols());
R = QR.matrixR();
printf("=[%d x %d]", R.rows(), R.cols());
for (i = 0; i < R.outerSize(); ++i) {
  for (SparseMatrix<double>::InnerIterator it(R, i); it; ++it) {
    printf("(%d %d %g)", it.row(), it.col(), it.value());
  }
}
------------------
Comment 1 Henrik Friberg 2014-06-28 03:05:58 UTC
Duh, I forgot to copy paste some lines:
-------------------
SparseMatrix<double> Q, R, AP;
long long int i;
-------------------

and
-------------------
AP = (A * QR.colsPermutation()).sparseView();
-------------------
Comment 2 Gael Guennebaud 2014-07-01 10:32:04 UTC
Indeed. TO be honest I thought the case columns>rows was explicitly not handled (doc+assert), but that was not the case. Anyway, this is now supported and covered by unit-tests:

https://bitbucket.org/eigen/eigen/commits/9c0c4bb28f32/
Changeset:   9c0c4bb28f32
User:        ggael
Date:        2014-07-01 10:24:46
Summary:     Fix bug 836: extend SparseQR to support more columns than rows.

and backported to 3.2: https://bitbucket.org/eigen/eigen/commits/657c9ed00ebc/
Comment 3 Henrik Friberg 2014-07-01 22:53:27 UTC
It means a lot to my time pressured ph.d. project -- thank you so much!!

Note You need to log in before you can comment on or make changes to this bug.