This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen
Bug 836 - Wrong results of SparseQR for non-square matrices
Summary: Wrong results of SparseQR for non-square matrices
Status: RESOLVED FIXED
Alias: None
Product: Eigen
Classification: Unclassified
Component: Sparse (show other bugs)
Version: 3.2
Hardware: All All
: Normal Wrong Result
Assignee: Nobody
URL:
Whiteboard:
Keywords:
Depends on:
Blocks:
 
Reported: 2014-06-28 03:02 UTC by Henrik Friberg
Modified: 2019-12-04 13:28 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!!
Comment 4 Nobody 2019-12-04 13:28:24 UTC
-- GitLab Migration Automatic Message --

This bug has been migrated to gitlab.com's GitLab instance and has been closed from further activity.

You can subscribe and participate further through the new bug through this link to our GitLab instance: https://gitlab.com/libeigen/eigen/issues/836.

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