C++ Eigen LU decomposition - sign of determinant

Question

I use C++ 14 and Eigen, I want to compute the determinant of a square matrix using LU decomposition. (link) but I have some issues: A is the main matrix (n size) and rez is the LU form of A.

PartialPivLU<MatrixXd> rez = PartialPivLU<MatrixXd>(A);
MatrixXd r1 = rez.matrixLU().triangularView<UpLoType::Upper>();
double det_r1 = 1;
for(int i=0;i<n;i++)
    det_r1 = det_r1 * r1(i,i);
cout<<det_r1<<endl;
cout<<A.determinant();

Determinant of r1 is the product of the elements from main diagonal. The problem is that det(A) isn't equal to det(r1).

For exemple det(A) = 500 and det(r1) = -500. The problem is about sign of r1, how can I get the sign?

I need to compute the determinant using the product of elements from main diagonal — dacian
As you have access to the indices of the permutation, you can calculate if this permutation is odd or even — Damien
@dacian That is what Eigen is doing internally. But it also keeps track on the number of permutations and thus correctly computes the sign as well (just check the source code to see what is happening). — chtz
Calculating the parity of a permutation needs a little effort, for example looking at the lengths of the cycles. You can find some information here: en.wikipedia.org/wiki/Parity_of_a_permutation — Damien

gspr gspr · Accepted Answer · 2020-07-20T09:49:52

The very documentation you link to says this is LU factorization with partial pivoting. That means you're decomposing A as PLU, where P is a permutation matrix. Its determinant is the sign of the permutation.

From the documentation you link to:

This class represents a LU decomposition of a square invertible matrix, with partial pivoting: the matrix A is decomposed as A = PLU where L is unit-lower-triangular, U is upper-triangular, and P is a permutation matrix.

In your case, the permutation matrix must represent an odd permutation, and therefore has determinant -1.

C++ Eigen LU decomposition - sign of determinant

1 Answers