FullPivLU cannot compute kernel of a null matrix.
Submitted by Matthieu Vigne
Assigned to Matthieu Vigne
Link to original bugzilla bug (#1617)
Version: 3.3 (current stable)
Description
Trying to compute the kernel of a nxn zero matrix using FullPivLU throws a floating point exception.
Minimal example:
#include <iostream>
#include <Eigen/Dense>
int main(int argc, char **argv)
{
// Compute kernel: this line will create a floating point exception.
Eigen::MatrixXd kernel = Eigen::Matrix2d::Zero().fullPivLu().kernel();
std::cout << "kernel" << kernel << std::endl;
return 0;
}
Expected output: the kernel of a null matrix is the entire space, the return kernel should be a basis of the space (typically: the identity 2x2 matrix)
Actual result: Floating point exception (core dumped)
This behavior is present at least since Eigen 3.2.8, and still present in the developement build (07fcdd14).
This exception comes from the triangular solving done inside LU/FullPivLu.h, l.685
m.topLeftCorner(rank(), rank())
.template triangularView<Upper>().solveInPlace(
m.topRightCorner(rank(), dimker)
);
For a zero matrix, rank() is zero, so we are solving for an empty matrix: this line is in fact has no effect, but TriangularSolver throws an exception.
To fix this, we can either:
- test before calling the TriangularSolver that rank is non-zero.
- inside the TriangularSolver, return immediately if trying to solve an empty problem.