Hello, i am not a very experienced user but i have to solve larger
Systems of linear equations with sparse, symmetric semi-positive definite matrices.

For this i wanted to use CHOLMOD. But i cant find a good tutorial on how to use the Cholmod-Support from within Eigen.

Even when i only include <unsupported/Eigen/SparseExtra/CholmodSupport.h>
i get tons of errors

like

::map has not been declared
::value has not been declared
cholmod_common does not name a type
...

I cant figure out what exactly i have to include and to write.

If someone could give me a minimum example how to
solve Ax = b using Cholmod via Eigen, id be really thankful.

Going back to sleep, eversleeping

First of all I recommend you to use Eigen 3.1.0-alpha2 which is much more mature regarding sparse stuff. Then have a look at the updated tutorial:

http://eigen.tuxfamily.org/dox-devel/Tu ... parse.html

For your first tries, you should simply use the built-in SimplicialLDLT solver which is quite fast for very sparse matrices. Then you can move to CholmodDecomposition.

I have to publish the resulting program. I am currently using cmake and some
FindEigen.cmake i have found in the web.

If my programm will use the new version, how can i check for it in cmake?

EDIT: i have found a FindEIgen in the downloads. I think this will be best option and it has a VERSION too.
Unfortunately this file is not installed when you install eigen via packagemanager.

So i have to copy this FindEigen3.cmake and add it to my cmake-folder ?

yes, you have to copy this file.

I am havin a little problem right now.
I am operating on matrices whcih should be (if calculated correctly)
positive-semi-definite.

The step i want to caluculate is

X = - [lambda * H1 + (1-lambda) H2 + Hc] ^-1 * (Grad1 - Grad2)

H stands for Hessian Matrices and Grads are gradients.

I am trying to solve the problem like this:

Code: Select all

Eigen::MatrixXd matrix=
         (  lambda     * enm1.calcHessian(current)
          + (1-lambda) * enm2.calcHessian(current)
          + scm.calcHessian(current)
          );

Eigen::VectorXd b = (enm1.calcGradient(current) - enm2.calcGradient(current));

Eigen::CholmodDecomposition<Eigen::SparseMatrix<double> > cd(matrix.sparseView());
cd.setMode(Eigen::CholmodLDLt);

X = - cd.solve(b);


But i get :

CHOLMOD warning: not positive definite

Do i use Eigen-Commands in some wrong way?
Or do i have to consider i calculate the matrices wrong or even the whole algorithm fails?

Thx

I see that you are using dense matrices. How large is 'matrix' ? Perhaps it would be better to simply use a dense solver:

LLDT<MatrixXd> ldlt(matrix);
X = ldlt.solve(-b);

Now, to check if 'matrix' is really positive definite, you could check its eigenvalues:

SelfAdjointEigenSolver<MatrixXd> eig(matrix,EigenvaluesOnly);
std::cout << eig.eigenvalue().minCoeff() << " " << eig.eigenvalue().maxCoeff() << std::endl;

I am getting confused more and more:

According to the paper I am reading, the usage of the sparse equation solver CHOLMOD is recommended.

Now i am calculating a simple matrix A (12*12):

1.7993 3.29773 2.19849 -2.85491 -3.29773 -2.19849 1.05561 0 0 0 0 0
3.29773 15.3975 6.91442 -3.29773 -5.60302 -3.29773 0 -9.79444 -3.61668 0 0 0
2.19849 6.91442 3.00486 -2.19849 -3.29773 -2.85491 0 -3.61668 -0.149955 0 0 0
-2.85491 -3.29773 -2.19849 5.15642 -1.13085 2.19849 -3.35712 4.42859 0 1.05561 0 0
-3.29773 -5.60302 -3.29773 -1.13085 22.4451 -0.318947 4.42859 -7.04761 0 0 -9.79444 3.61668
-2.19849 -3.29773 -2.85491 2.19849 -0.318947 3.40959 0 0 -0.404729 0 3.61668 -0.149955
1.05561 0 0 -3.35712 4.42859 0 5.15642 -1.13085 -2.19849 -2.85491 -3.29773 2.19849
0 -9.79444 -3.61668 4.42859 -7.04761 0 -1.13085 22.4451 0.318947 -3.29773 -5.60302 3.29773
0 -3.61668 -0.149955 0 0 -0.404729 -2.19849 0.318947 3.40959 2.19849 3.29773 -2.85491
0 0 0 1.05561 0 0 -2.85491 -3.29773 2.19849 1.7993 3.29773 -2.19849
0 0 0 0 -9.79444 3.61668 -3.29773 -5.60302 3.29773 3.29773 15.3975 -6.91442
0 0 0 0 3.61668 -0.149955 2.19849 3.29773 -2.85491 -2.19849 -6.91442 3.00486

The real problem has matrices about 1000*1000.
But nevertheless:

Vector b:

-1.31284
4.36438
0.798374
2.12229
7.0887
-0.798374
-2.12229
-7.0887
-0.798374
1.31284
-4.36438
0.798374

The Cholmodsupport-Module never gave any good solution. Now i have like adviced calculated the eigenvalues.

-1.60069
-1.286
-0.450127
-6.69925e-16
-7.43596e-17
5.43958e-15
0.225471
7.023
7.81154
22.1422
32.9951
35.5649

Now i have tried some of the offered solvers:

USING LU().solve()
-1.33253
1.00044
0.931639
-0.388127
1.25973
-0.024051
-1.20242
0.763206
-0.024051
-0.258022
1.0225
0.931639

USING FullPivLu()
-0.202744
-0.026851
0.47121
0.741654
0.232443
-0.48448
-0.0726381
-0.264082
-0.48448
0.87176
-0.00478803
0.47121

USING ldlt() should not be working (not semi-definite)
-1.0745
-0.0220629
-1.02478e-15
-0.130106
0.237231
-0.95569
-0.944398
-0.259294
-0.95569
0
1.72654e-15
8.94704e-17

FullPivHouseHolderQr
0
-0.0220629
0.95569
0.944398
0.237231
0
0.130106
-0.259294
-2.98802e-16
1.0745
0
0.95569

I wonder what the best solution may be.

Given the matrix to maple i get the solution:

-0.537246
-0.011033
0.4952
..

yet another solution. But it may be that some digits behind the comma are not print so that Maple has less precision and therefor another solution.

So how should i continue from here on?

Your problem is underdetermined, or in other words not full rank, so there no a single solution but there exist a space of solutions.You need to add some additional constraints. For instance you can use a SVD to get the minimal norm solution. For a 1000x1000 matrix, it is still ok to use dense algorithms:

x = JacobiSVD<MatrixXd>(matrix, ComputeThinU|ComputeThinV).solve(b);

Moreover, it seems your matrix has a band structure. If that's really the case, the best for you would be use an Eigen::BandMatrix to assemble the matrix, and then call the lapack DSBEV routine. (we currently don't have optimized algorithm for band matrices).

The Matrix A

1.7993 3.29773 2.19849 -2.85491 -3.29773 -2.19849 1.05561 0 0 0 0 0
3.29773 15.3975 6.91442 -3.29773 -5.60302 -3.29773 0 -9.79444 -3.61668 0 0 0
2.19849 6.91442 3.00486 -2.19849 -3.29773 -2.85491 0 -3.61668 -0.149955 0 0 0
-2.85491 -3.29773 -2.19849 5.15642 -1.13085 2.19849 -3.35712 4.42859 0 1.05561 0 0
-3.29773 -5.60302 -3.29773 -1.13085 22.4451 -0.318947 4.42859 -7.04761 0 0 -9.79444 3.61668
-2.19849 -3.29773 -2.85491 2.19849 -0.318947 3.40959 0 0 -0.404729 0 3.61668 -0.149955
1.05561 0 0 -3.35712 4.42859 0 5.15642 -1.13085 -2.19849 -2.85491 -3.29773 2.19849
0 -9.79444 -3.61668 4.42859 -7.04761 0 -1.13085 22.4451 0.318947 -3.29773 -5.60302 3.29773
0 -3.61668 -0.149955 0 0 -0.404729 -2.19849 0.318947 3.40959 2.19849 3.29773 -2.85491
0 0 0 1.05561 0 0 -2.85491 -3.29773 2.19849 1.7993 3.29773 -2.19849
0 0 0 0 -9.79444 3.61668 -3.29773 -5.60302 3.29773 3.29773 15.3975 -6.91442
0 0 0 0 3.61668 -0.149955 2.19849 3.29773 -2.85491 -2.19849 -6.91442 3.00486

Has as far as i can see full rank.
Even if some eigenvalues are suspiciously small.
But nevertheless the system of equations should have an unique solution.

If i am wrong (it is my first time that i do such stuff), plz tell me where
i err.
Thx

The banded structure like in this example is always occuring in some form but not a single band but more, something like:

X 0 0 Y Z 0
0 X 0 0 Y Z
...

well, with double 1e-15 compared to 35 means 0. I tested with the printed values and I got:

singular vals: 35.565 32.9951 22.1422 7.81155 7.023 1.60069 1.286 0.450117 0.225467 3.57001e-05 4.49997e-06 7.00281e-07
eigen vals: -1.60069 -1.286 -0.450117 -4.49997e-06 -7.00281e-07 3.57001e-05 0.225467 7.023 7.81155 22.1422 32.9951 35.565

FullPivLU -0.537247 -0.0110334 0.495247 0.407148 0.248263 -0.460445 -0.407148 -0.248263 -0.460445 0.537247 0.0110334 0.495247
PartialPivLU -0.537247 -0.0110334 0.495247 0.407148 0.248263 -0.460445 -0.407148 -0.248263 -0.460445 0.537247 0.0110334 0.495247
JacobiSVD -0.537247 -0.0110334 0.495247 0.407148 0.248263 -0.460445 -0.407148 -0.248263 -0.460445 0.537247 0.0110334 0.495247
LDLT -0.537247 -0.0110334 0.495247 0.407148 0.248263 -0.460445 -0.407148 -0.248263 -0.460445 0.537247 0.0110334 0.495247

resp errs: 1.96104e-15 4.14665e-15 1.31497e-09 2.25897e-15

The clamped version is indeed positive definite and easy to solve with any solver. Perhaps you could try to regularize your problem by adding epsilon to the diagonal.... Our built-in sparse LDLT solver can do that for you:

SimplicialLDLT<SparseMatrix<double> > ldlt;
ldlt.setShift(epsilon);

ldlt.compute(A);
if(ldlt.info()!=Success)
increase epsilon and try again !

x = ldtd.solve(b);