Hi,

I'm having trouble with eigen sparse matrix llt with the Cholmod backend.

I'm running eigen 2.0.12.

basically i get the wrong answer.

here is some code showing my problem, Sorry if this is the wrong forum for this question.

#include <iostream>

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

int main( int argc, char * argv[])
{
// USE EIGEN Sparse Matrix
typedef Eigen::SparseMatrix<double> SparseMatrixType;
typedef Eigen::DynamicSparseMatrix<double> DynamicSparseMatrixType;
typedef Eigen::VectorXd DenseVectorType;
typedef Eigen::MatrixXd DenseMatrixType;

DynamicSparseMatrixType m_tmp = DynamicSparseMatrixType(8,8);

for (unsigned int i=0;i<8;++i)
{
m_tmp.coeffRef(i,i) = 1.0;
}
m_tmp.coeffRef(0,1) = -0.25;
m_tmp.coeffRef(0,3) = -0.25;
m_tmp.coeffRef(1,0) = -0.25;
m_tmp.coeffRef(1,2) = -0.25;
m_tmp.coeffRef(2,1) = -0.25;
m_tmp.coeffRef(2,4) = -0.25;
m_tmp.coeffRef(2,1) = -0.25;
m_tmp.coeffRef(2,4) = -0.25;
m_tmp.coeffRef(3,0) = -0.25;
m_tmp.coeffRef(3,5) = -0.25;
m_tmp.coeffRef(4,2) = -0.25;
m_tmp.coeffRef(4,7) = -0.25;
m_tmp.coeffRef(5,3) = -0.25;
m_tmp.coeffRef(5,6) = -0.25;
m_tmp.coeffRef(6,5) = -0.25;
m_tmp.coeffRef(6,7) = -0.25;
m_tmp.coeffRef(7,4) = -0.25;
m_tmp.coeffRef(7,6) = -0.25;
SparseMatrixType m_sparse(m_tmp);

DenseVectorType rhs = DenseVectorType::Zero(8);
rhs[1] = 0.25;
rhs[3] = 0.25;
rhs[4] = 0.25;
rhs[6] = 0.25;

DenseMatrixType m_dense = m_sparse.toDense();
DenseVectorType x_dense(8);
DenseVectorType x_sparse = rhs;

Eigen::SparseLLT<SparseMatrixType,Eigen::Cholmod> sparseLLT(m_sparse);
Eigen::LLT<DenseMatrixType> denseLLT(m_dense);

sparseLLT.solveInPlace(x_sparse);
denseLLT.solve(rhs, &x_dense);

std::cout << "Sparse Matrix : " << std::endl << m_sparse << std::endl << std::endl;
std::cout << "Dense Matrix : " << std::endl << m_dense << std::endl << std::endl;
std::cout << "rhs : " << std::endl << rhs << std::endl << std::endl;
std::cout << "x_dense : " << std::endl << x_dense << std::endl << std::endl;
std::cout << "x_sparse : " << std::endl << x_sparse << std::endl << std::endl;

return EXIT_SUCCESS;
}

The output xdense and x_Sparse output is...

x_dense :
0.1667
0.3333
0.1667
0.3333
0.3333
0.1667
0.3333
0.1667

x_sparse :
0.4444
0.5556
0.4444
0.5556
0.5556
0.4444
0.5556
0.4444

I appreciate any help / advice.
Luke

There are lots of problems in what you're doing here

First of all, your matrix:

Code: Select all

for (unsigned int i=0;i<8;++i)
{
m_tmp.coeffRef(i,i) = 1.0;
}
m_tmp.coeffRef(0,1) = -0.25;
m_tmp.coeffRef(0,3) = -0.25;
m_tmp.coeffRef(1,0) = -0.25;
m_tmp.coeffRef(1,2) = -0.25;
m_tmp.coeffRef(2,1) = -0.25;
m_tmp.coeffRef(2,4) = -0.25;
m_tmp.coeffRef(2,1) = -0.25;
m_tmp.coeffRef(2,4) = -0.25;
m_tmp.coeffRef(3,0) = -0.25;
m_tmp.coeffRef(3,5) = -0.25;
m_tmp.coeffRef(4,2) = -0.25;
m_tmp.coeffRef(4,7) = -0.25;
m_tmp.coeffRef(5,3) = -0.25;
m_tmp.coeffRef(5,6) = -0.25;
m_tmp.coeffRef(6,5) = -0.25;
m_tmp.coeffRef(6,7) = -0.25;
m_tmp.coeffRef(7,4) = -0.25;
m_tmp.coeffRef(7,6) = -0.25;


is certainly NOT positive definite, since it has zeros on the diagonal. So you can't use a LLT decomposition with that matrix. Actually, since it has zero diagonal, it can't be even positive semidefinite.

But there's another point... this matrix is singular. Indeed we have this relation between columns:

Code: Select all

col(1) - col(3) - col(4) + col(6) = 0

so no matter what decomposition you use, the solution of solve() is not unique --- there exist different equally valid solutions.

Thanks for the response.

The matrix is full rank, it has ones down the diagonal?

Printing it yields...

1 -0.25 0 -0.25 0 0 0 0
-0.25 1 -0.25 0 0 0 0 0
0 -0.25 1 0 -0.25 0 0 0
-0.25 0 0 1 0 -0.25 0 0
0 0 -0.25 0 1 0 0 -0.25
0 0 0 -0.25 0 1 -0.25 0
0 0 0 0 0 -0.25 1 -0.25
0 0 0 0 -0.25 0 -0.25 1

matlab returns the eigen values as
1.5000
1.3536
1.3536
1.0000
1.0000
0.6464

Are you running ubuntu? I've been having all types of trouble getting eigen 2.0.12 to recognise the different backends like superLU etc.

thanks,
Luke

My bad. I didn't see that at the beginning you filled the diagonal with 1's.

Sorry. I don't know much about the sparse module, I hope someone can help you.

np, I appreciate any help i can get thanks.

Luke

you must fill the lower triangular part only