Different results when naming a matrix-vector product used in triangular solve, vs passing using it directly
Submitted by Sebastian Sylvan
Assigned to Nobody
Link to original bugzilla bug (#646)
Version: 3.2
Description
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()).triangularViewEigen::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