This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen
Bug 1479 - LDLT thinks isPositive() is true for an indefinite matrix
Summary: LDLT thinks isPositive() is true for an indefinite matrix
Status: RESOLVED FIXED
Alias: None
Product: Eigen
Classification: Unclassified
Component: Cholesky (show other bugs)
Version: 3.3 (current stable)
Hardware: All All
: Normal Unknown
Assignee: Nobody
URL:
Whiteboard:
Keywords:
Depends on:
Blocks:
 
Reported: 2017-10-21 03:11 UTC by Hongkai Dai
Modified: 2019-12-04 17:14 UTC (History)
2 users (show)



Attachments
Test code to reproduce the output. (640 bytes, text/x-c++src)
2017-10-21 03:11 UTC, Hongkai Dai
no flags Details

Description Hongkai Dai 2017-10-21 03:11:55 UTC
Created attachment 798 [details]
Test code to reproduce the output.

I tested a matrix with both positive and negative eigen values, and LDLT thinks the matirx is positive definite, when I call "isPositive()" function. The test code is
#include <Eigen/Core>
#include <Eigen/Cholesky>
#include <Eigen/Eigenvalues>
#include <iostream>

int main() {
	Eigen::Matrix4d M;
	M << 1, 2, 0, 1,
	     2, 4, 0, 2,
	     0, 0, 0, 1,
	     1, 2, 1, 1;
	Eigen::LDLT<Eigen::Matrix4d> ldlt(M);
	std::cout << "ldlt.info() == Eigen::Success? " << (ldlt.info() == Eigen::Success) << "\n\n";
	std::cout << "ldlt.isPositive()? " << ldlt.isPositive() << "\n\n";

	Eigen::SelfAdjointEigenSolver<Eigen::Matrix4d> es(M);
	std::cout << "Eigen values:\n" << es.eigenvalues().transpose() << "\n\n";
	Eigen::Vector4d x;
	x << 0, 0, -1, 1;
	std::cout << "x'*M*x = " << x.dot(M * x) << "\n\n";
	return 0;
}

The output is

ldlt.info() == Eigen::Success? 1

ldlt.isPositive()? 1

Eigen values:
  -0.924984 1.07385e-16    0.896688      6.0283

x'*M*x = -1


Attached is the code.

I checked the `D` vector in LDLT, it is all non-negative. Shouldn't LDLT.isPositive() return false for an in-definite matrix?
Comment 1 Gael Guennebaud 2017-11-16 17:05:10 UTC
OK, the true problem is that it returned success whereas the factorization failed. With this fix, it "properly" return NumericalIssue.

https://bitbucket.org/eigen/eigen/commits/2ff106034f40/ (devel)
https://bitbucket.org/eigen/eigen/commits/67719139abc3/ (3.3)

Actually, with a proper permutation it is possible to find a LDL^T factorization of this matrix, for instance, the following works fine (if we disable permutation in Eigen::LDLT):

1 1 2 1
1 0 0 0
2 0 4 2
1 0 2 1

Unfortunately, I don't think it is possible to quickly find the proper permutation automatically during the factorization. Currently, it simply permutes the matrix according to the magnitude of the diagonal entries. I also tried a right-looking variant:

      Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
      Block<MatrixType,Dynamic,Dynamic> A22(mat,k+1,k+1,rs,rs);
      RealScalar realAkk = numext::real(mat.coeffRef(k,k));
      bool pivot_is_valid = (abs(realAkk) > RealScalar(0));
      if((rs>0) && pivot_is_valid)
        A21 /= realAkk;
      if(rs>0)
        A22.template selfadjointView<Lower>().rankUpdate(A21,-realAkk);

but that does not help. For instance, if we start with the largest pivot '4', then the remaining matrix block A22 become:

0 1 0
1 0 0
0 0 0

and we are stuck.
Comment 2 Nobody 2019-12-04 17:14:58 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/1479.

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