Hello,
i am trying to solve some equation in my program, but somehow i fail.

Matrix A, it is singular:

1.27689 -2.5395 -1.83667 -1.24699 2.66815 1.9483 -0.0299056 -0.128647 -0.11163 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
-2.5395 6.26238 4.64894 2.66815 -5.70897 -4.16873 -0.128647 -0.553409 -0.480206 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
-1.83667 4.64894 3.46073 1.9483 -4.16873 -3.04404 -0.11163 -0.480206 -0.416685 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
-1.24699 2.66815 1.9483 5.56727 0.934014 1.85454 -4.21979 -3.37617 -3.60454 -0.100495 -0.225992 -0.198303 1.17097e-16 1.12024e-16 1.44979e-16 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
2.66815 -5.70897 -4.16873 0.934014 8.91838 7.4986 -3.37617 -2.70121 -2.88393 -0.225992 -0.508204 -0.445938 1.12024e-16 1.36683e-16 1.58206e-16 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1.9483 -4.16873 -3.04404 1.85454 7.4986 6.51435 -3.60454 -2.88393 -3.079 -0.198303 -0.445938 -0.391301 1.44979e-16 1.58206e-16 2.19187e-16 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
-0.0299056 -0.128647 -0.11163 -4.21979 -3.37617 -3.60454 4.53091 3.45347 3.83012 -0.0367389 0.482967 0.364393 -0.148682 -0.204563 -0.291084 -0.0957949 -0.227054 -0.187256 0 0 0 0 0 0 0 0 0 0 0 0
-0.128647 -0.553409 -0.480206 -3.37617 -2.70121 -2.88393 3.45347 10.4233 8.99874 0.482967 -6.34904 -4.79028 -0.204563 -0.281446 -0.400485 -0.227054 -0.538166 -0.443835 0 0 0 0 0 0 0 0 0 0 0 0
-0.11163 -0.480206 -0.416685 -3.60454 -2.88393 -3.079 3.83012 8.99874 8.04582 0.364393 -4.79028 -3.61422 -0.291084 -0.400485 -0.569872 -0.187256 -0.443835 -0.366039 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 -0.100495 -0.225992 -0.198303 -0.0367389 0.482967 0.364393 5.55455 0.86017 5.05918 -5.15524 -0.786519 -4.93531 -0.262074 -0.330626 -0.289961 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 -0.225992 -0.508204 -0.445938 0.482967 -6.34904 -4.79028 0.86017 7.39435 6.355 -0.786519 -0.119997 -0.752965 -0.330626 -0.41711 -0.365808 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 -0.198303 -0.445938 -0.391301 0.364393 -4.79028 -3.61422 5.05918 6.355 9.0511 -4.93531 -0.752965 -4.72477 -0.289961 -0.365808 -0.320816 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 1.17097e-16 1.12024e-16 1.44979e-16 -0.148682 -0.204563 -0.291084 -5.15524 -0.786519 -4.93531 5.63967 2.52686 5.77283 -0.165586 -1.23503 -0.321136 -0.170168 -0.300747 -0.225305 0 0 0 0 0 0 0 0 0
0 0 0 1.12024e-16 1.36683e-16 1.58206e-16 -0.204563 -0.281446 -0.400485 -0.786519 -0.119997 -0.752965 2.52686 10.1446 3.94686 -1.23503 -9.21161 -2.39522 -0.300747 -0.531526 -0.398192 0 0 0 0 0 0 0 0 0
0 0 0 1.44979e-16 1.58206e-16 2.19187e-16 -0.291084 -0.400485 -0.569872 -4.93531 -0.752965 -4.72477 5.77283 3.94686 6.21575 -0.321136 -2.39522 -0.622808 -0.225305 -0.398192 -0.298305 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 -0.0957949 -0.227054 -0.187256 -0.262074 -0.330626 -0.289961 -0.165586 -1.23503 -0.321136 4.34029 4.00159 5.64092 -3.65651 -2.03942 -4.363 -0.124943 -0.146972 -0.296195 -0.0353833 -0.0224902 -0.183372 0 0 0
0 0 0 0 0 0 -0.227054 -0.538166 -0.443835 -0.330626 -0.41711 -0.365808 -1.23503 -9.21161 -2.39522 4.00159 11.4915 6.1033 -2.03942 -1.13749 -2.43347 -0.146972 -0.172885 -0.348418 -0.0224902 -0.0142951 -0.116555 0 0 0
0 0 0 0 0 0 -0.187256 -0.443835 -0.366039 -0.289961 -0.365808 -0.320816 -0.321136 -2.39522 -0.622808 5.64092 6.1033 8.16816 -4.363 -2.43347 -5.20601 -0.296195 -0.348418 -0.702172 -0.183372 -0.116555 -0.950322 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 -0.170168 -0.300747 -0.225305 -3.65651 -2.03942 -4.363 3.95415 2.58149 5.392 -0.0529296 -0.334591 -0.64385 -0.00565959 -0.00156088 0.0750008 -0.0688874 0.0948219 -0.234842
0 0 0 0 0 0 0 0 0 0 0 0 -0.300747 -0.531526 -0.398192 -2.03942 -1.13749 -2.43347 2.58149 3.91505 6.55777 -0.334591 -2.11509 -4.07005 -0.00156088 -0.000430481 0.0206848 0.0948219 -0.13052 0.323254
0 0 0 0 0 0 0 0 0 0 0 0 -0.225305 -0.398192 -0.298305 -4.363 -2.43347 -5.20601 5.392 6.55777 15.1308 -0.64385 -4.07005 -7.83198 0.0750008 0.0206848 -0.99391 -0.234842 0.323254 -0.800593
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.124943 -0.146972 -0.296195 -0.0529296 -0.334591 -0.64385 0.673131 1.26525 -0.598172 -0.410552 -0.999672 1.71395 -0.084707 0.215987 -0.175731
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.146972 -0.172885 -0.348418 -0.334591 -2.11509 -4.07005 1.26525 5.27285 -0.202986 -0.999672 -2.43415 4.17338 0.215987 -0.550725 0.448081
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.296195 -0.348418 -0.702172 -0.64385 -4.07005 -7.83198 -0.598172 -0.202986 16.054 1.71395 4.17338 -7.1553 -0.175731 0.448081 -0.364568
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.0353833 -0.0224902 -0.183372 -0.00565959 -0.00156088 0.0750008 -0.410552 -0.999672 1.71395 5.07459 -3.8709 -0.656636 -4.623 4.89463 -0.94894
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.0224902 -0.0142951 -0.116555 -0.00156088 -0.000430481 0.0206848 -0.999672 -2.43415 4.17338 -3.8709 7.63109 -5.0822 4.89463 -5.18222 1.0047
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.183372 -0.116555 -0.950322 0.0750008 0.0206848 -0.99391 1.71395 4.17338 -7.1553 -0.656636 -5.0822 9.29431 -0.94894 1.0047 -0.194784
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.0688874 0.0948219 -0.234842 -0.084707 0.215987 -0.175731 -4.623 4.89463 -0.94894 4.77659 -5.20544 1.35951
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0948219 -0.13052 0.323254 0.215987 -0.550725 0.448081 4.89463 -5.18222 1.0047 -5.20544 5.86346 -1.77603
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.234842 0.323254 -0.800593 -0.175731 0.448081 -0.364568 -0.94894 1.0047 -0.194784 1.35951 -1.77603 1.35995


Vector b:

0.0418035
0.0181618
0.0231441
0.543176
0.616413
0.662352
-0.253417
0.101262
0.0824684
-0.117774
-0.40587
-0.354135
-0.209276
-0.13494
-0.357029
0.487509
0.130048
0.660695
-0.412528
0.174613
0.260247
-0.211814
-1.01704
-0.0273045
0.427267
0.10854
-0.88716
-0.294945
0.408815
-0.0632786

Solution by A.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(b):

-1.34437e+13
1.38791e+13
-9.37773e+12
-1.05878e+13
1.69186e+13
-1.17124e+13
-1.20033e+13
2.03772e+12
3.88285e+12
-7.58083e+12
5.30146e+12
2.97824e+09
-1.18956e+13
-1.74857e+11
5.38271e+12
-2.70619e+12
-6.17293e+11
2.34597e+12
3.10091e+12
-3.88272e+12
-9.9441e+11
1.09513e+13
2.72233e+12
-5.07223e+12
3.62175e+12
1.34083e+13
-5.95255e+11
-7.72361e+12
3.83486e+12
5.29692e+12

Those values are much too big. Moreover calculating Ax does not result in b anymore (could be something numerical due to the big values).

Both Maple and Lapack give me better solutions.

One solution by Maple:
(1)=-3079.063356,(2)=1657.893814,(3)=475.7515633,(5)=514.9918500,(4)=-1532.578568,(7)=-149.4807078,(6)=3030.744242,(10)=1510.166901,(11)=-2910.859715,(8)=-1354.858660,(9)=3162.828063,(15)=-1353.613507,(14)=-569.4767024,(13)=7611.207394,(12)=5392.446937,(21)=998.6261607,(20)=531.1276857,(23)=1213.982458,(22)=-495.3072496,(16)=5011.242259,(17)=-411.1325992,(18)=-622.0137422,(19)=2551.804690,(30)=3231.142703,(29)=2258.659255,(28)=-2680.652207,(26)=1676.944447,(27)=621.4374982,(24)=894.1760006,(25)=-2760.810802])"

Dont know why the String export messes up the order.

So, how can i solve this and similar equations and get some "good" solution?

What about x = A.colPivHouseholderQr().solve(b);

I think that the problem is that the very small singular values are not clamped to 0 in JacobiSVD.

x = A.colPivHouseholderQr().solve(b);

Unfortunately not working well. The first 3 entries of the solution are zero now. But that should not be the case.

hm, here I get with JacobiSvd:
-3079.06 1657.89 475.752 -1532.58 514.992 3030.74 -149.481 -1354.86 3162.83 1510.17 -2910.86 5392.45 7611.21 -569.477 -1353.61 5011.24 -411.133 -622.014 2551.8 531.128 998.626 -495.307 1213.98 894.176 -2760.81 1676.94 621.437 -2680.65 2258.66 3231.14

and with colPivHouseholderQr:

-3079.06 1657.89 475.752 -1532.58 514.992 3030.74 -149.481 -1354.86 3162.83 1510.17 -2910.86 5392.45 7611.21 -569.477 -1353.61 5011.24 -411.133 -622.014 2551.8 531.128 998.626 -495.307 1213.98 894.176 -2760.81 1676.94 621.437 -2680.65 2258.66 3231.14

Make sure you are using double and not float, float are bad for singular problems.

Hi,
I also have a problem of calling Eigen to solve Ax=b, where the size of A is 113 x 113.
A is a symmetric sparse complex value-based matrix. I set x=1.0 and compute the b.
However, after calling both the dense and sparse solvers (SuperLU, Pardiso, etc.) from Eigen. I cannot obtain a solution of x=1.0 anymore. I tried the matlab x=A\b, it works.

The condition number of A is about 1.0E+8.

I think Eigen's direct solver might be not accurate enough.


Zhengyong

no need to use a sparse solver for such a small matrix. Please can you output the matrix to a file and send it to send me for testing?

ggael wrote:no need to use a sparse solver for such a small matrix. Please can you output the matrix to a file and send it to send me for testing?

Dear Guennebaud,

I send the matrix named: Krr_0.dat to your email already. Thanks for your help.

Best wishes
Zhengyong

btw, do you use the compiler flag

Code: Select all

-O3

?

manuels wrote:btw, do you use the compiler flag

Code: Select all

-O3

?

Yes, I used -O3 and the matrix I sent to you is in matlab format with starting index=1.

Dear all,

I tried to input the matrix again and using Eigen, now it works.


Thanks a lot for your help.


Best wishes
Zhengyong

As for my problem, i will post something more next week. By now i dont have a clou what is going wrong.

The Matrix and vector are from

std::cout << Eigen::MatrixXd

maybe the precision is not good enough when printed and internally it calculates something more complex?

renzh: your matrix works fine here with all Eigen's solvers:

MatrixXd A;
VectorXd x, b;
SparseMatrix<double> S;
loadMarket(S,"Krr_0.dat");
int n = S.cols();
A = MatrixXd(S);
x.setOnes(n);
b = A*x;

x.setZero();
x = SimplicialLDLT<SparseMatrix<double> >(S).solve(b);
//x = A.ldlt().solve(b);
//x = A.partialPivLu().solve(b);
//x = A.colPivHouseholderQr().solve(b);
//x = A.jacobiSvd(ComputeFullU|ComputeFullV).solve(b);

Eversleeping: you might output your matrix with full precision:

std::cout << std::setprecision (numeric_limits<double>::digits10 + 1) << mat;

ggael wrote:Eversleeping: you might output your matrix with full precision:

std::cout << std::setprecision (numeric_limits<double>::digits10 + 1) << mat;

Yes, the previous matrix "Krr_0.dat" was not outputted with a full precision. I also get Eigen working on the non full-precision outputted Krr_0.dat from Eigen. Now, I tried the " std::setprecision (numeric_limits<double>::digits10 + 1)" and std::setprecision(16) (which was already resent to you) to output again. The condition number becomes 1e+17. It is closely singular. And matlab also cannot work.

The background of the matrix is from the a sub-domain in the domain decomposition. It might be that the geometry of this sub-domain is not good enough.


Thanks for your help.

Best wishes
Zhengyong

ggael wrote:renzh: your matrix works fine here with all Eigen's solvers:

MatrixXd A;
VectorXd x, b;
SparseMatrix<double> S;
loadMarket(S,"Krr_0.dat");
int n = S.cols();
A = MatrixXd(S);
x.setOnes(n);
b = A*x;

x.setZero();
x = SimplicialLDLT<SparseMatrix<double> >(S).solve(b);
//x = A.ldlt().solve(b);
//x = A.partialPivLu().solve(b);
//x = A.colPivHouseholderQr().solve(b);
//x = A.jacobiSvd(ComputeFullU|ComputeFullV).solve(b);

Hi, my matrix is a complex one, could you please try the second matrix I sent you by email with MatrixEC format.

Best wishes
Zhengyong