[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Bug-gsl] LU_inverse bug
|
From: |
Georges Szegezdi |
|
Subject: |
[Bug-gsl] LU_inverse bug |
|
Date: |
Mon, 11 Oct 2004 16:05:38 +0200 |
|
User-agent: |
Mozilla Thunderbird 0.8 (Windows/20040913) |
Hi, I'm using the GSL library to do the pseudo-inverse calcul of a matrix (and
then do some other stuff on it like transpose etc..).
I have a problem with some type of matrices, I reproduce the example above.
When inverting a matrix, the result sometimes comes out to NaN or Inf. I tried
both of the example above in Matlab and Octave and it worked, giving a warning tho.
I think LU_invert should catch these errors(or impossible calculations) and
reply with a very high number instead of a NaN, so we can continue the
calculation with the inverted matrix. Or return an error code or something, I've
looked around but haven't found anything about such return code(apologize me if
that exists).
Here is the code that I use:
gsl_matrix *Z = gsl_matrix_calloc (N,N);
gsl_blas_dgemm(CblasTrans,CblasNoTrans,1.0,anotherMatrix,anotherMatrix,0.0,Z);
gsl_matrix *inverse = gsl_matrix_calloc (N,N);
gsl_permutation *p = gsl_permutation_calloc(N);
int s;
gsl_linalg_LU_decomp(Z, p, &s);
gsl_linalg_LU_invert(Z, p, inverse);
The values and result of the matrices are:
Z:
1 0 0 0 0 0 0 0.988097 -0.821691
1.07345e-05 0.000741653
0 0 0 0 0 0 0 0 0 0 0
0 0 3 0 0 0 0 2.97242 -1.80443
0.00179485 0.000916459
0 0 0 2 0 0 0 1.97459 -1.69657
0.00118993 0.00113353
0 0 0 0 1 0 0 0.981981 -0.803057
0.000930131 0.000391623
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 2 1.91503 -1.38932 1.00122
1.00186
0.988097 0 2.97242 1.97459 0.981981 0 1.91503 8.67102
-6.39055 0.932407 0.99256
-0.821691 0 -1.80443 -1.69657 -0.803057 0
-1.38932 -6.39055 4.81479 -0.748694 -0.648082
1.07345e-05 0 0.00179485 0.00118993 0.000930131 0
1.00122 0.932407 -0.748694 1.00199 0.00115457
0.000741653 0 0.000916459 0.00113353 0.000391623 0
1.00186 0.99256 -0.648082 0.00115457 1.00186
inverse:
NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
NaN
NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
NaN
NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
NaN
NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
NaN
NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
NaN
NaN NaN NaN NaN NaN Inf NaN NaN NaN NaN
NaN
1.87121e+08 0 1.77284e+08 1.91665e+08 1.89854e+08 0
8.18867e+09 -1.40127e+08 5.18902e+07 -8.00467e+09 -8.00764e+09
-3.21066e+06 -0 -3.04114e+06 -3.28855e+06 -3.25729e+06 -0
-1.40127e+08 2.40206e+06 -893389 1.3697e+08 1.37021e+08
1.21369e+06 -0 1.14531e+06 1.2432e+06 1.23044e+06 -0
5.18902e+07 -893389 356305 -5.06987e+07 -5.072e+07
-1.82888e+08 0 -1.73278e+08 -1.87329e+08 -1.8556e+08 0
-8.00467e+09 1.3697e+08 -5.06987e+07 7.82482e+09 7.82773e+09
-1.82958e+08 0 -1.73344e+08 -1.87401e+08 -1.85631e+08 0
-8.00764e+09 1.37021e+08 -5.072e+07 7.82773e+09 7.83063e+09
If I add tiny little noise (about 0.0001) to the 0 and 1's, then Z and inverse
are the following, and the result can be computed or used later on.
Z:
1.00196 0.000371977 0.00153351 0.00237999 0.00160393 0.000778553
0.00173078 0.145228 0
.181509 0.00139875 0.001183
0.000371977 3.28104e-06 0.00235466 0.00119582 7.42628e-05
2.57241e-06 0.000684671 -1.93966e
-05 0.000254357 8.50162e-05 0.000603617
0.00153351 0.00235466 3.00054 0.00249093 0.00166448 0.00132093
0.00275591 0.0597847 -
0.0701054 0.00276492 0.00303812
0.00237999 0.00119582 0.00249093 2.00327 0.00152111 0.00140153
0.00159401 -0.404662 0
.72661 0.00201412 0.00115836
0.00160393 7.42628e-05 0.00166448 0.00152111 1.00194
0.000226316 0.000708809 0.085846 0
.101533 0.000619092 0.000642688
0.000778553 2.57241e-06 0.00132093 0.00140153 0.000226316
3.46488e-06 0.00125373 0.0001779
1 0.000336786 0.000780943 0.000476505
0.00173078 0.000684671 0.00275591 0.00159401 0.000708809
0.00125373 2.00203 0.448053 -
0.537788 1.00138 1.00199
0.145228 -1.93966e-05 0.0597847 -0.404662 0.085846
0.00017791 0.448053 0.220735-
0.215039 0.289155 0.159107
0.181509 0.000254357 -0.0701054 0.72661 0.101533 0.000336786
-0.537788 -0.215039 0
.517096 -0.129895 -0.407587
0.00139875 8.50162e-05 0.00276492 0.00201412 0.000619092
0.000780943 1.00138 0.289155 -
0.129895 1.00097 0.00103852
0.001183 0.000603617 0.00303812 0.00115836 0.000642688
0.000476505 1.00199 0.159107 -
0.407587 0.00103852 1.00168
inverse:
-2.11121e+14 1.37696e+17 -1.35357e+14 4.36086e+14 -1.15437e+14
-3.6991e+16 -5.1383e+14 1.76733e+
15 -3.70722e+14 -2.75925e+13 1.72982e+13
1.37981e+17 -8.99931e+19 8.84663e+16 -2.85002e+17 7.54453e+16
2.41745e+19 3.37537e+17 -1.15504e
+18 2.42288e+17 1.63131e+16 -1.30239e+16
-1.3788e+14 8.99291e+16 -8.84179e+13 2.84732e+14 -7.53876e+13
-2.41454e+16 -3.50761e+14 1.15403e+
15 -2.42091e+14 -2.79554e+12 2.65045e+13
4.26692e+14 -2.78286e+17 2.73498e+14 -8.81619e+14 2.33318e+14
7.48099e+16 9.81945e+14 -3.57263e
+15 7.49341e+14 1.12455e+14 2.16622e+13
-1.15066e+14 7.50469e+16 -7.37697e+13 2.37687e+14 -6.29163e+13
-2.01628e+16 -2.77811e+14 9.63263e+
14 -2.02055e+14 -1.72821e+13 7.18693e+12
-3.52257e+16 2.29732e+19 -2.25715e+16 7.28093e+16 -1.92627e+16
-6.18101e+18 -7.51066e+16 2.95014e+
17 -6.18707e+16 -1.5257e+16 -7.75471e+15
-2.59225e+15 1.69241e+18 -1.67716e+15 5.29793e+15 -1.41516e+15
-4.43568e+17 -1.88197e+16 2.15435e+
16 -4.53382e+15 1.22034e+16 1.27402e+16
1.74131e+15 -1.13568e+18 1.11622e+15 -3.59751e+15 9.52145e+14
3.05232e+17 4.08137e+15 -1.45788e
+16 3.05791e+15 3.84631e+14 1.41962e+13
-3.67718e+14 2.39828e+17 -2.35734e+14 7.5963e+14 -2.01065e+14
-6.44441e+16 -8.76872e+14 3.07846e+
15 -6.45728e+14 -6.6191e+13 1.20177e+13
2.05606e+15 -1.34271e+18 1.33343e+15 -4.19025e+15 1.12197e+15
3.49594e+17 1.75463e+16 -1.70546e
+16 3.59225e+15 -1.23051e+16 -1.27279e+16
2.09856e+15 -1.37043e+18 1.36067e+15 -4.27812e+15 1.14522e+15
3.57054e+17 1.76354e+16 -1.74106e
+16 3.66691e+15 -1.22851e+16 -1.2717e+16
Here is a last example, I also had NaN on that matrix, I don't know why tho.
this one doesn't use pure 0 or 1's.
Z:
297929 3805.55 -3825.03 3825.03 -3825.03 -3825.03 1241.64
3805.55 50.4217 -50.7088 50.7088 -50.7088 -50.7088 14.0719
-3825.03 -50.7088 51 -51 51 51 -13.7362
3825.03 50.7088 -51 51 -51 -51 13.7362
-3825.03 -50.7088 51 -51 51 51 -13.7362
-3825.03 -50.7088 51 -51 51 51 -13.7362
1241.64 14.0719 -13.7362 13.7362 -13.7362 -13.7362 102.012
inverse:
NaN NaN NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN NaN NaN
NaN NaN -Inf NaN Inf NaN NaN
5.81559e+12 -5.24055e+16 2.06566e+30 -0 -0 -2.06566e+30
2.82898e+14
0 0 -2.81475e+14 0 0 2.81475e+14 0
I don't know if it is a real bug or a math issue with matrices, please let me
know.
I have looked into the mailing forum but haven't found anything about that
topic.
Thanks.
G
- [Bug-gsl] LU_inverse bug,
Georges Szegezdi <=