Re: [eigen] Cholesky Decompositon

[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]


Dear Tommaso,

in the basic use of the LLT, the input matrix is copied to an internal memory (in advanced use, the input matrix can be used as internal memory to avoid this copy). The decomposition is then done in place, from the left to the right, i.e. the column of the internal memory are sequentially replaced by the columns of the matrix L. The algorithm works only on the lower triangular part of the matrix and thus the strictly upper triangular part of the internal memory is left untouched.
If a non-positive pivot is encountered, the decomposition stops, leaving the internal memory in its current state, meaning that the rightmost part of the internal memory will be left untouched and thus still contains the elements of the original matrix.
The matrixL() function returns a lower triangular view on the internal memory.
I can't tell what's happening in case other problems (which ones?) arise.

The code below can help illustrate this behavior.

Best regards,
Adrien


#include <Eigen/Core>
#include <Eigen/Cholesky>
#include <Eigen/QR>

#include <iostream>

using namespace Eigen;

int main()
{
  int n = 5;
  //random eigen values, all positive
  VectorXd s = VectorXd::Random(n).cwiseAbs();
  s[0] = -s[0]; //make one eigen value negative. Change which one to see different behavior

  //Generate a random orthogonal matrix
  MatrixXd R = MatrixXd::Random(n, n);
  HouseholderQR<MatrixXd> qr(R);
  const auto& Q = qr.householderQ();

  //M is a symmetric matrix with all eigen values positive but one
  MatrixXd M = Q * s.asDiagonal().toDenseMatrix() * Q.transpose();
  std::cout << "M=\n" << M << "\n\n";

  auto llt = M.llt(); //perform Cholesky decomposition
  std::cout << "info: " << llt.info() << "\n";

  std::cout << "L=\n" << llt.matrixL().toDenseMatrix() << "\n\n";
  std::cout << "Internal memory:\n" << llt.matrixLLT() << "\n\n";

  std::cout << "What elements are different between L and the lower part of M (1 when different, 0 when equal):\n";
  std::cout << (M.template triangularView<Lower>().toDenseMatrix().cwiseEqual(llt.matrixL().toDenseMatrix())).select(MatrixXd::Zero(n, n), MatrixXd::Ones(n, n)) << "\n\n";
  std::cout << "What elements are different between the internal memory and M (1 when different, 0 when equal):\n";
  std::cout << M.cwiseEqual(llt.matrixLLT()).select(MatrixXd::Zero(n, n), MatrixXd::Ones(n, n)) << std::endl;
}

On Tue, Apr 14, 2020 at 10:36 PM Tommaso Ferrari <tomferri93@xxxxxxxxx> wrote:
Dear Adrien, 
thank you for your kind and precise reply. 
I tried to invoke the info method and I get exactly what you reported. I would therefore like to ask you the following. 

In the llt() method, if a problem of any nature were encountered, an exception would not be raised, but the matrix L would be terminated with the elements of the initial matrix.
Is this correct?

Thank you in advance

Il giorno mar 14 apr 2020 alle ore 15:00 Adrien Escande <adrien.escande@xxxxxxxxx> ha scritto:
Dear Tommaso,

calling the method info() on your LLT object in this case should return Eigen::ComputationInfo::NumericalIssue. It's not because the code didn't crash that you can say the decomposition arrived at a result. Success is conveyed by info() and there is no exception raised upon failure.
From a quick look a the code, llt() stops prematurely if it encounters a pivot that is <= 0. The part of the resulting L matrix up to this point will be "correct", but all the rows after should simply be (part of) the rows of your initial matrix.

Best regards,
Adrien

On Tue, Apr 14, 2020 at 7:35 PM Tommaso Ferrari <tomferri93@xxxxxxxxx> wrote:
Dear all,
I was analyzing Cholesky's decomposition algorithm on a non positive definite matrix. It is a 3*3 matrix, whose eigenvalues are -29.5, 2, 30.5. 
I noticed that the llt() method produces a result even if, obviously, the reconstruction of the starting matrix is not correct (due to the fact that the input matrix is not positive definite). 
I would like to know if, in general, given a non positive definite matrix, the llt() method still arrives at a result or if there may be exceptions or errors. 
How does the decomposition proceed, having as input a non positive definite matrix? Are pseudo method (pseudo determinant, pseudo inverse) used?
Thanks for your attention and response

Tommaso 


Mail converted by MHonArc 2.6.19+ http://listengine.tuxfamily.org/