Hello,
I am trying to iteratively solve a Ax=b problem . I am using a gauss-newton algorithm, and the normal equation matrix A is of the form:
A=|U W Ct 0 |
|Wt V Bt Gt|
|C B D 0 |
|0 G 0 0 |
Sub matrices U, V, D are block diagonal matrices and Ct is for C.transpose() and so on... A is about 11000*11000 and is quite sparse (80% of zeros)
I've used 3 different algorithm to solve Ax=b with success (2 dense and a sparse one): colPivHouseholderQr, LU, and SparseLU. They all give me the same results with same rms. Good.
Next, i've tried SimplicialLDLT. In the documentation A is required to be symetric and positive definite. A is effectively symmetric but i don't know if it is positive definite. Anyway, it works great, info.() return Success, same results as above and same rms, and it's FAST!
Next, i've tried the dense version : LDLT. The results are differents and not good, but with same rms, info() return success() and both isPositive() and isNegative() return false!

Here is my question:
I'm planning to use SimplicialLDLT because it's the fastest and accurate. But, could i trust the results, as the dense version give me wrong results ?

Thanks.

hm, that's rather odd that the sparse version does work while the dense does not because the dense version has the advantage of doing numerical pivoting. So the dense version should be much more accurate. Can you save the SparseMatrix and right-hand-side as follows:

Code: Select all

#include <unsupported/Eigen/SparseExtra>
...
SparseMatrix<double> A;
VectorXd b;
...
Eigen::saveMarket(A, "filename.mtx",Eigen::Symmetric);
Eigen::saveMarketVector(B, "filename_b.mtx");


so that we can have a look. thanks.

Hi, and thanks for you answer.

Here is the link to both matrix A and vector b:
https://www.transfernow.net/fr/0gaoacxc3vf6
The matrix is smaller than in my previous post( 2000x2000 if i am correct) but same behaviour.

I don't know if it is usefull for you, but here are some additionals informations:
- I use sparseView() to get the sparse matrix A.
- I use Eigen 3.2.0
- It's a constrained non-linear least square problem. In my first post the G matrix is a constraint matrix, and the last 7 values in b are Lagrangian multipliers (to avoid the rank deficiency).

Thanks.

Alright, the dense LDLT is fixed in the devel branch: https://bitbucket.org/eigen/eigen/commits/f240bbe84f94/.

In your case I'd recommend to use SimplicialLDLT, check the residual, e.g.:

(A*x-b).norm()/b.norm()<NumTraits<double>::dummy_precision(),

and if it fails, fallback to SparseLU. The reason is that your matrix is not positive definite, and since SimplicialLDLT does not perform numerical pivoting it might fail in some cases.

Thank you very much!