This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen
Bug 1617 - FullPivLU cannot compute kernel of a null matrix.
Summary: FullPivLU cannot compute kernel of a null matrix.
Status: RESOLVED FIXED
Alias: None
Product: Eigen
Classification: Unclassified
Component: LU (show other bugs)
Version: 3.3 (current stable)
Hardware: All All
: Normal Crash
Assignee: Matthieu Vigne
URL:
Whiteboard:
Keywords: test-needed
Depends on:
Blocks:
 
Reported: 2018-10-30 14:05 UTC by Matthieu Vigne
Modified: 2019-12-04 18:04 UTC (History)
3 users (show)



Attachments
Minimal patch for the current bug. (2.74 KB, patch)
2018-10-31 19:53 UTC, Matthieu Vigne
no flags Details | Diff

Description Matthieu Vigne 2018-10-30 14:05:55 UTC
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 (9584fac85b3f).

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.
Comment 1 Christoph Hertzberg 2018-10-30 14:41:30 UTC
Good catch!

I'd prefer having TriangularSolver (in fact any solver) return immediately (after asserting that the sizes are compatible) for 0x0-matrices.
Comment 2 Matthieu Vigne 2018-10-31 19:53:16 UTC
Created attachment 889 [details]
Minimal patch for the current bug.

I agree fixing the solver is more generic. However from what I could see, there are separate implementation of triangular solvers in both Core/SolveTriangular and SparseCore/TriangularSolver.

I've made a minimal patch for the current bug: it only modifies Core/SolveTriangular solveInPlace function. This is enough to get FullPivLU.kernel() to work (unit test added).

If you see a better, more generic way to patch this, feel free to offer. This is my first time patching Eigen and I didn't go to far into working out the internals, but I can look deeper.
Comment 3 Gael Guennebaud 2018-11-01 13:45:56 UTC
Thank you for the patch. Looks ok for dense triangular solve.

https://bitbucket.org/eigen/eigen/commits/359f78f1be47/
https://bitbucket.org/eigen/eigen/commits/7d79b21a2ed8/ (3.3)

Ultimately, we should likely apply similar fixes for sparse triangular solves and all solve functions.
Comment 4 Gael Guennebaud 2019-01-16 14:25:41 UTC
https://bitbucket.org/eigen/eigen/commits/99e9c99cc552/
Summary:     Bug 1617: add unit tests for empty triangular solve.
Comment 5 Nobody 2019-12-04 18:04:37 UTC
-- 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/1617.

Note You need to log in before you can comment on or make changes to this bug.