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 <[email protected]>
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 <
> [email protected]> 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 <[email protected]>
>> 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
>>>
>>
