CODEWITHSUNDEEP
c++eigen
user1559897
user1559897

Reputation: 1484

What is the correct way to solve a least square equation with eigen?

I have the following simple example to perform a least square, but got the following assertion error.

Assertion failed: (i>=0) && ( ((BlockRows==1) && (BlockCols==XprType::ColsAtCompileTime) && i

What is the correct way to do this?

typedef Eigen::ArrayXd arr;
typedef Eigen::ArrayXXd arr2;

arr2 A(3, 3);
arr B(3);
A << 1, 0, 0, 0, 1, 0, 0, 0, 1;
B << 1, 2, 3;
auto x = A.matrix().colPivHouseholderQr().solve(B.matrix());

Upvotes: 1

Views: 1245

Answers (1)

ggael
ggael

Reputation: 29205

As said in my comment, the problem is that x is an abstract expression storing a reference to the QR object, but right after the last line it will reference a dead object and then anything can happen!

More precisely, A.matrix().colPivHouseholderQr() creates a temporary object, let's call it tmp_qr. Then tmp_qr.solve(B) creates another object which becomes x of type Solve<...>. This object essentially stores two references: one to tmp_qr, and one to B. Right after this line, the temporary object tmp_qr is deleted, so the Solve<...> object x has one dead reference. It is like a pointer referencing a deleted buffer. Finally, if you x later, like:

VectorXd y = x;

operator= will trigger the solve operation using the QR decomposition referenced by x and the right-hand-side B referenced by x too, but, wait... the QR decomposition object has been deleted, so at best you'll get a segfault.

So the solution is to write:

VectorXd x = A.matrix().colPivHouseholderQr().solve(B.matrix());

See Eigen's doc for more details on how auto is dangerous if you don't know what you get.

Upvotes: 3

Related Questions