While trying to use QR factorization to solve linear least squares I've encountered a bug. In the final step I get wildly different results depending on whether or not I name an intermediate vector. Just introducing a temporary fixes it. There's no aliasing anywhere. Here's a code snippet that demonstrates the issue (using VS2012 and VS2013 preview): // A is a NxM linear matrix (I used 100x3 for plane fitting) // y is Nx1 matrix representing measurements. // We're minimizing the square of r, where r = Ab-y auto QR = A.fullPivHouseholderQr(); // Solve for R1 b' = Q1^t y', where R1 and Q1^t hold the first 3 rows. // where b' = Pb, P from the QR factorization // This gives least squares solution. auto R1 = QR.matrixQR().topRows(A.cols()).triangularView<Eigen::Upper>(); auto Q1t = QR.matrixQ().adjoint().topRows(A.cols()); // BUG shows up in the following. // Let's solve for b three different different ways. auto tmp = Q1t*y; auto bprime = R1.solve(tmp); Eigen::Vector3f b = QR.colsPermutation().inverse() * bprime; // THIS returns the right result (notice the 'tmp' above) auto bprime2 = R1.solve(Q1t*y); Eigen::Vector3f b2 = QR.colsPermutation().inverse() * bprime2; // THIS returns wildly incorrect results Eigen::Vector3f bprime3 = R1.solve(Q1t*y); // Force it to eval to vec3 Eigen::Vector3f b3 = QR.colsPermutation().inverse() * bprime3; // THIS returns the correct results again
(In reply to comment #0) > auto bprime2 = R1.solve(Q1t*y); > Eigen::Vector3f b2 = QR.colsPermutation().inverse() * bprime2; // THIS returns > wildly incorrect results I assume this is related to bug 505. bprime2 is an expression which depends on a product and the product is evaluated into a temporary before solve is called. But solve() itself is only evaluated when bprime2 is used one line later and at that time the product temporary is not valid anymore. I assume your first way only works by accident. If you really like to write auto instead of the actual type for your temporaries, you should add an (...).eval(); to the corresponding expressions if they involve products. You can, btw, directly write Eigen::Vector3f b4 = QR.solve(y);
Okay, so I should evaluate products for now. Kinda of a big problem for readability (a lot of expressions involve products). Does eval() actually evaluate the whole thing, or does it doe some kind of shallow evaluation? I.e. will it still do the expression template goodness (e.g. if I only end up accessing a single element from the matrix)? Re: auto, it seems like a godsend for stuff like this. Now you can build complicated expression with named temporaries that still get the expression template goodness (merely naming a temporary doesn't force evaluation). Ah, thanks for the hint on QR.solve(y), I wasn't sure if it would produce a least-squares minimal solution or not (as an aside, perhaps docs could be updated to mention this, they could also be updated to describe how to get the R matrix, since this is not obvious either).
BTW: Eigen::Vector3f b4 = QR.solve(y); Does not return the correct result (or anything close to it).
(In reply to comment #2) > Okay, so I should evaluate products for now. Kinda of a big problem for > readability (a lot of expressions involve products). The easiest solution would be to write temporaries into actual Eigen::Vector## objects, no need to write eval() then. > Does eval() actually evaluate the whole thing, or does it doe some kind of > shallow evaluation? I.e. will it still do the expression template goodness > (e.g. if I only end up accessing a single element from the matrix)? It evaluates everything and generally, it's not a good idea to access a single element from complex expressions. > Re: auto, it seems like a godsend for stuff like this. Now you can build > complicated expression with named temporaries that still get the expression > template goodness (merely naming a temporary doesn't force evaluation). What's the point of a temporary, if it's not evaluated? If you don't want to evaluate it, don't make a temporary for it but directly insert the expression into the next expression. If you need the temporary more than once, you should evaluate it before you use it (otherwise, it gets evaluated twice, unless your compiler is really smart). > Ah, thanks for the hint on QR.solve(y), I wasn't sure if it would produce a > least-squares minimal solution or not (as an aside, perhaps docs could be > updated to mention this, they could also be updated to describe how to get the > R matrix, since this is not obvious either). Actually, we should just add a matrixR() method (e.g. ColPivQR has one as well). (In reply to comment #3) > Eigen::Vector3f b4 = QR.solve(y); > > Does not return the correct result (or anything close to it). Hm, sorry. I was assuming and expecting that it gives the least squares solution. I can look into that, but I'm afraid not too soon. Maybe Gael can explain if the current result is intended.
> If you don't want to evaluate it, don't make a temporary for it but directly insert the expression into the next expression. But then you end up with *massive* expressions which are hard to read. It's much better to be able to name the individual parts of an expression as you go, for readability (without losing the template expression goodness). > Hm, sorry. I was assuming and expecting that it gives the least squares solution. It may well be supposed to. The value returned isn't just wrong in some subtle way, it seems to be returning garbage. I get this back: {-107374176., -107374176., -107374176.}
The testcases for FullPivHouseholderQR indeed do not seem to cover injective but non-surjective matrices (as they appear in LS problems), but that's actually an unrelated issue. For the readability thing, I don't really agree. If the expression gets very huge, it is almost certainly worthwhile to explicitly evaluate some intermediate results into actual matrices -- that's what Eigen does as well at the moment. (I don't deny that the current behavior with `auto` temporaries still is a bug).
*** This bug has been marked as a duplicate of bug 505 ***
(In reply to comment #3) > BTW: > Eigen::Vector3f b4 = QR.solve(y); > > Does not return the correct result (or anything close to it). Works for me: #include <Eigen/Dense> #include <iostream> using namespace Eigen; typedef Matrix<double, Dynamic, 2> MatrixX2d; typedef Matrix<double, Dynamic, 3> MatrixX3d; int main() { int n = 100; MatrixX3d D(n,3); VectorXd f(n); Vector3d x, u; // Initial plane equation u << 11.11, -0.201005, 1.8; // Generate random samples MatrixX2d pts = MatrixX2d::Random(n, 2); // Evaluate the inital model, and add some noise D.col(0).fill(1); D.rightCols<2>() = pts; f = D * u + 0.05*VectorXd::Random(n); std::cout << "Solution: " << u.transpose() << "\n"; // Solve x = D.colPivHouseholderQr().solve(f); std::cout << "Found (ColQR): " << x.transpose() << " (" << (D*x-f).norm() << ")\n"; x = D.fullPivHouseholderQr().solve(f); std::cout << "Found (FullQR): " << x.transpose() << " (" << (D*x-f).norm() << ")\n"; x = MatrixXd(D).jacobiSvd(ComputeThinU|ComputeThinV).solve(f); std::cout << "Found (SVD): " << x.transpose() << " (" << (D*x-f).norm() << ")\n"; x = (D.adjoint()*D).ldlt().solve(D.adjoint()*f); std::cout << "Found (LDLt): " << x.transpose() << " (" << (D*x-f).norm() << ")\n"; }
Odd, I've pasted a complete test case below, using JacobiSVD as the standard, and seeing how the QR decomposition based version fares. The assert consistently fires (the QR version always returns {-107374176., -107374176., -107374176.} #include <Eigen/Eigen> #include <cassert> static const int NUM_POINTS = 100; int main(int argc, char* argv[]) { auto y = Eigen::VectorXf::Random(NUM_POINTS); auto A = Eigen::MatrixXf::Random(NUM_POINTS, 3); Eigen::Vector3f b = A.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(y); auto QR = A.fullPivHouseholderQr(); Eigen::Vector3f b2 = QR.solve(y); assert(b.isApprox(b2, 0.1f)); return 0; }
(In reply to comment #9) > auto y = Eigen::VectorXf::Random(NUM_POINTS); > auto A = Eigen::MatrixXf::Random(NUM_POINTS, 3); It's again basically the same problem as above. A is just an expression, which gets evaluated into a temporary, and gets invalid after the next `;` Even, if all expression would be evaluated as late as possible, you would get different results here.
No, the same problem happens if you get rid of "auto".
Just to be clear, this still fires the assert (and again, the QR result returns what appears to be garbage). #include <Eigen/Eigen> #include <cassert> static const int NUM_POINTS = 100; int main(int argc, char* argv[]) { Eigen::VectorXf y = Eigen::VectorXf::Random(NUM_POINTS); Eigen::MatrixXf A = Eigen::MatrixXf::Random(NUM_POINTS, 3); Eigen::Vector3f b = A.jacobiSvd(Eigen::ComputeThinU|Eigen::ComputeThinV).solve(y); Eigen::Vector3f b2 = A.fullPivHouseholderQr().solve(y); assert(b.isApprox(b2, 0.1f)); return 0; }
alright, I've never really looked at the solve method of full-piv QR, and indeed it was not doing a LS solving. Fixed: https://bitbucket.org/eigen/eigen/commits/6814a13155ec/ Changeset: 6814a13155ec User: ggael Date: 2013-08-20 11:52:48 Summary: Make FullPivHouseholderQR::solve returns the least-square solution instead of aborting if no exact solution exist
I can confirm that this fixes my issue. I wasn't even aware that the solve method could fail (and I'm still not quite sure how to check!). Thanks for the quick fix! I guess I can live with naming all my return values for now.
-- 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/646.