search

adamaulia / efficient-java-matrix-library

Automatically exported from code.google.com/p/efficient-java-matrix-library
0 stars 0 forks source link

Erroneous calculation of the determinant of a matrix #38

Closed GoogleCodeExporter closed 9 years ago

GoogleCodeExporter commented 9 years ago 
The calculation of the determinant of a Matrix seems to increase a lot when the 
dimension of a square matrix increases. Way(!) beyond any appropriate tolerance.

I would aspect at most a tolerance about 1e-15 for the following test, while 
the error value seems to increase to infinity with increasing dimension.

The following shows the results for the difference between the determinant of a 
random matrix and its transpose (det(A) - det(A^T)). The correct value would be 
0 for each given case.

1. column: Dimension
2. column: Error using Colt
3. column: Error using EJML
4. column: Error using EJML for det(A) - det((A^t)^t), to differentiate whether 
det() or ^t seems to be erroneous.

1   0.0         0.0         0.0
2   -5.551115123125783E-17  0.0         0.0
3   2.7755575615628914E-17  1.3877787807814457E-17  0.0
4   2.7755575615628914E-17  1.3877787807814457E-17  0.0
5   -6.938893903907228E-18  -6.938893903907228E-17  0.0
6   -6.938893903907228E-18  3.642919299551295E-17   0.0
7   -6.938893903907228E-18  -3.469446951953614E-18  0.0
8   -2.7755575615628914E-17 3.469446951953614E-18   0.0
9   -2.7755575615628914E-17 1.9081958235744878E-17  0.0
10  2.7755575615628914E-17  -5.551115123125783E-17  0.0
11  2.7755575615628914E-17  3.122502256758253E-17   0.0
12  6.938893903907228E-18   -6.288372600415926E-18  0.0
13  6.938893903907228E-18   0.0         0.0
14  6.938893903907228E-18   1.610040226140974E-17   0.0
15  -3.122502256758253E-17  9.367506770274758E-17   0.0
16  -1.3877787807814457E-17 -2.7755575615628914E-17 0.0
17  -1.3877787807814457E-17 -1.0625181290357943E-17 0.0
18  1.3877787807814457E-17  -7.806255641895632E-17  0.0
19  -1.3877787807814457E-17 5.551115123125783E-17   0.0
20  0.0         -1.734723475976807E-18  0.0
21  0.0         -6.938893903907228E-17  0.0
22  6.938893903907228E-18   0.0         0.0
23  0.0         -9.43689570931383E-16   0.0
24  -1.734723475976807E-18  5.551115123125783E-17   0.0
25  6.938893903907228E-17   1.8318679906315083E-15  0.0
26  1.734723475976807E-18   2.7755575615628914E-17  0.0
27  2.7755575615628914E-17  -2.7755575615628914E-15 0.0
28  2.7755575615628914E-17  1.3877787807814457E-15  0.0
29  9.71445146547012E-17    -3.552713678800501E-15  0.0
30  6.938893903907228E-18   6.217248937900877E-15   0.0
31  3.469446951953614E-18   -1.5987211554602254E-13 0.0
32  5.551115123125783E-17   -3.552713678800501E-14  0.0
33  5.551115123125783E-17   4.618527782440651E-14   0.0
34  -1.3877787807814457E-17 -5.790923296444817E-13  0.0
35  2.7755575615628914E-17  1.7053025658242404E-13  0.0
36  5.551115123125783E-17   2.2737367544323206E-13  0.0
37  5.551115123125783E-17   -4.831690603168681E-13  0.0
38  -5.204170427930421E-18  2.3447910280083306E-13  0.0
39  -5.551115123125783E-17  1.0800249583553523E-12  0.0
40  -5.551115123125783E-17  1.8189894035458565E-12  0.0
41  8.326672684688674E-17   1.1084466677857563E-12  0.0
42  8.326672684688674E-17   1.5916157281026244E-12  0.0
43  -1.1102230246251565E-16 -9.094947017729282E-12  0.0
44  -2.7755575615628914E-17 -7.457856554538012E-11  0.0
45  -2.7755575615628914E-17 2.5579538487363607E-13  0.0
46  8.326672684688674E-17   1.6825651982799172E-11  0.0
47  -6.938893903907228E-18  2.1373125491663814E-10  0.0
48  -9.540979117872439E-18  -3.166496753692627E-8   0.0
49  2.7755575615628914E-17  2.6193447411060333E-10  0.0
50  -6.938893903907228E-18  -3.725290298461914E-9   0.0
51  2.7755575615628914E-17  -1.0739313438534737E-8  0.0
52  1.0408340855860843E-17  4.470348358154297E-7    0.0
53  1.0408340855860843E-17  1.0710209608078003E-8   0.0
54  -5.551115123125783E-17  9.575160220265388E-9    0.0
55  -2.7755575615628914E-17 -8.032657206058502E-9   0.0
56  0.0         4.842877388000488E-8    0.0
57  -1.0408340855860843E-17 5.9604644775390625E-8   0.0
58  -1.1102230246251565E-16 0.0         0.0
59  1.3877787807814457E-17  1.430511474609375E-6    0.0
60  6.938893903907228E-18   -3.4332275390625E-5 0.0
61  1.3877787807814457E-17  0.0         0.0
62  2.7755575615628914E-17  -3.039836883544922E-6   0.0
63  -1.3877787807814457E-17 7.32421875E-4       0.0
64  -1.3877787807814457E-17 2.09808349609375E-5 0.0
65  6.938893903907228E-18   -0.00140380859375   0.0
66  -8.673617379884035E-19  -0.001613616943359375   0.0
67  2.7755575615628914E-17  -6.103515625E-5     0.0
68  -6.938893903907228E-18  0.0010986328125     0.0
69  5.551115123125783E-17   -0.0028076171875    0.0
70  1.1492543028346347E-17  -0.029296875        0.0
71  4.163336342344337E-17   0.0078125       0.0
72  3.469446951953614E-18   0.25            0.0
73  -4.163336342344337E-17  6.46875         0.0
74  0.0         -2.4375         0.0
75  -5.551115123125783E-17  -0.5625         0.0
76  1.3877787807814457E-17  2.4375          0.0
77  0.0         1.5         0.0
78  -1.3877787807814457E-17 1120.0          0.0
79  3.8163916471489756E-17  296.0           0.0
80  -1.3877787807814457E-17 -4160.0         0.0
81  5.551115123125783E-17   1088.0          0.0
82  5.095750210681871E-18   -52736.0        0.0
83  0.0         34176.0         0.0
84  4.163336342344337E-17   5632.0          0.0
85  6.938893903907228E-18   24832.0         0.0
86  5.551115123125783E-17   293072.0        0.0
87  -5.551115123125783E-17  229376.0        0.0
88  5.551115123125783E-17   1638400.0       0.0
89  1.3877787807814457E-17  -851968.0       0.0
90  8.348356728138384E-18   1245184.0       0.0
91  -1.3877787807814457E-17 851968.0        0.0
92  0.0         -7.5497472E7        0.0
93  -2.7755575615628914E-17 -6.3963136E7        0.0
94  2.0816681711721685E-17  4.9283072E7     0.0
95  -6.938893903907228E-18  8.9915392E7     0.0
96  5.551115123125783E-17   7.38197504E8        0.0
97  -4.163336342344337E-17  -9.663676416E9      0.0
98  0.0         -2.952790016E9      0.0
99  -1.1102230246251565E-16 -2.10453397504E11   0.0

Code snippet to reproduce:

for(int i = 1; i < 100; i++) {
  // Dimension
  System.out.print(i + "\t");

  // Error using Colt
  DoubleFactory2D doubleFactory2D = DoubleFactory2D.dense;
  DoubleMatrix2D matrixCern = doubleFactory2D.random(4,4);
  Algebra algebra = new Algebra();
  System.out.print(algebra.det(matrixCern) - algebra.det(algebra.transpose(matrixCern)) + "\t");

  // Error using EJML ...
  Random random = new Random();
  DenseMatrix64F matrix = RandomMatrices.createRandom(i, i, random);
  System.out.print(CommonOps.det(matrix) - CommonOps.det(CommonOps.transpose(matrix, null)) + "\t");
  // ... for det(A) - det((A^t)^t), to differentiate whether det() or ^t seems to be erroneous.
  System.out.println(CommonOps.det(matrix) - CommonOps.det(CommonOps.transpose(CommonOps.transpose(matrix, null), null)));
}

Original issue reported on code.google.com by os...@united-forum.de on 2 Oct 2013 at 9:37

GoogleCodeExporter commented 9 years ago 
Thanks for the bug report.  I see a few problems with the test, the largest one 
being that the Colt matrix is always 4x4 while EJML keeps on growing.

I also did some tests in Octave, which is a wrapper around LAPACK, and it 
appears to have a similar response to EJML.

octave:15> a=rand(3,3);
octave:16> det(a)-det(a')
ans = -2.7756e-17

octave:11> a=rand(60,60);
octave:12> det(a)-det(a')
ans =  5.7220e-06

a=rand(80,80);
octave:14> det(a)-det(a')
ans =  96

Could you make the following change? 

DoubleMatrix2D matrixCern = doubleFactory2D.random(4,4);

to

DoubleMatrix2D matrixCern = doubleFactory2D.random(i,i);

Original comment by peter.ab...@gmail.com on 2 Oct 2013 at 10:30

GoogleCodeExporter commented 9 years ago 
Thanks for the fast response and sorry for this silly error on my side.

The results between Colt and EJML are quite similar now. It was actually my 
first concern to check whether this is a problem unique to EJML or more some 
kind of problem that also affects other libraries.

1   0.0         0.0         0.0
2   6.938893903907228E-18   0.0         0.0
3   0.0         -6.938893903907228E-18  0.0
4   6.071532165918825E-18   -1.734723475976807E-18  0.0
5   2.7755575615628914E-17  5.551115123125783E-17   0.0
6   6.938893903907228E-18   -2.2985086056692694E-17 0.0
7   3.0357660829594124E-18  -2.2551405187698492E-17 0.0
8   1.3877787807814457E-17  6.938893903907228E-17   0.0
9   -2.7755575615628914E-17 0.0         0.0
10  4.3368086899420177E-19  -6.938893903907228E-18  0.0
11  6.938893903907228E-17   2.0816681711721685E-17  0.0
12  -1.3877787807814457E-17 0.0         0.0
13  1.6479873021779667E-17  5.204170427930421E-17   0.0
14  -2.7755575615628914E-17 -1.942890293094024E-16  0.0
15  -1.3877787807814457E-17 -4.163336342344337E-17  0.0
16  -3.469446951953614E-17  -3.469446951953614E-17  0.0
17  -5.551115123125783E-17  4.5102810375396984E-17  0.0
18  2.986976985197565E-17   5.692061405548898E-18   0.0
19  2.0816681711721685E-17  7.979727989493313E-17   0.0
20  2.2551405187698492E-17  -1.214306433183765E-17  0.0
21  -1.3877787807814457E-17 1.249000902703301E-16   0.0
22  4.5102810375396984E-17  -2.0816681711721685E-17 0.0
23  6.661338147750939E-16   9.992007221626409E-16   0.0
24  -4.884981308350689E-15  -2.4424906541753444E-15 0.0
25  7.216449660063518E-16   2.220446049250313E-16   0.0
26  -3.2751579226442118E-15 -1.6514567491299204E-15 0.0
27  1.942890293094024E-15   -6.106226635438361E-16  0.0
28  9.71445146547012E-16    0.0         0.0
29  -1.7208456881689926E-15 2.9976021664879227E-15  0.0
30  4.6629367034256575E-15  4.440892098500626E-16   0.0
31  9.103828801926284E-15   -2.7977620220553945E-14 0.0
32  -2.6645352591003757E-14 4.618527782440651E-14   0.0
33  2.5951463200613034E-15  9.506284648352903E-15   0.0
34  -1.0231815394945443E-12 -4.405364961712621E-13  0.0
35  3.268496584496461E-13   2.5579538487363607E-13  0.0
36  -2.4868995751603507E-14 1.1723955140041653E-13  0.0
37  9.379164112033322E-13   -5.115907697472721E-13  0.0
38  2.2737367544323206E-13  5.684341886080801E-13   0.0
39  5.542233338928781E-13   -4.973799150320701E-13  0.0
40  -1.8971491044794675E-12 -1.3358203432289883E-12 0.0
41  -6.423306331271306E-12  -2.4442670110147446E-12 0.0
42  3.865352482534945E-12   4.774847184307873E-12   0.0
43  6.548361852765083E-11   -3.637978807091713E-12  0.0
44  -2.9558577807620168E-12 -1.1937117960769683E-11 0.0
45  -9.276845958083868E-11  2.5283952709287405E-10  0.0
46  -2.3283064365386963E-10 6.402842700481415E-10   0.0
47  -1.0986695997416973E-9  -9.64064383879304E-10   0.0
48  4.0745362639427185E-9   -8.149072527885437E-10  0.0
49  3.92901711165905E-10    -3.055902197957039E-10  0.0
50  -1.3969838619232178E-9  3.259629011154175E-9    0.0
51  4.6566128730773926E-9   1.3504177331924438E-8   0.0
52  9.313225746154785E-9    2.3283064365386963E-9   0.0
53  1.30385160446167E-8 -2.60770320892334E-8    0.0
54  6.48200511932373E-7 4.470348358154297E-8    0.0
55  -6.332993507385254E-8   4.0978193283081055E-8   0.0
56  5.245208740234375E-6    4.5299530029296875E-6   0.0
57  -1.5832483768463135E-7  -1.2293457984924316E-7  0.0
58  -5.066394805908203E-7   -8.642673492431641E-7   0.0
59  -2.682209014892578E-6   4.351139068603516E-6    0.0
60  2.384185791015625E-6    -5.900859832763672E-6   0.0
61  -2.8520822525024414E-5  -1.621246337890625E-5   0.0
62  1.41143798828125E-4 2.09808349609375E-4 0.0
63  3.814697265625E-6   -8.20159912109375E-5    0.0
64  8.678436279296875E-5    8.58306884765625E-6 0.0
65  4.2724609375E-4     0.004180908203125   0.0
66  -5.6743621826171875E-5  1.1968612670898438E-4   0.0
67  7.62939453125E-6    -1.02996826171875E-4    0.0
68  -0.0142822265625    -0.0169677734375    0.0
69  0.3828125       0.4921875       0.0
70  -0.03515625     -0.09375        0.0
71  0.052490234375      -0.006103515625     0.0
72  0.01953125      -0.0556640625       0.0
73  0.0015606880187988281   0.0026330947875976562   0.0
74  -28.59375       -7.75           0.0
75  1.09375         0.046875        0.0
76  -12.0           -22.9375        0.0
77  -6.0            -62.0           0.0
78  64.0            832.0           0.0
79  -2272.0         640.0           0.0
80  148.0           -88.0           0.0
81  -64.0           -816.0          0.0
82  -13056.0        -20480.0        0.0
83  -1408.0         -2848.0         0.0
84  -6784.0         14848.0         0.0
85  -139264.0       700416.0        0.0
86  -58368.0        86016.0         0.0
87  786432.0        -32768.0        0.0
88  311296.0        638976.0        0.0
89  1015808.0       -81920.0        0.0
90  655360.0        1.2189696E7     0.0
91  -6.5011712E7        -1.6777216E7        0.0
92  1.92937984E8        2.3986176E8     0.0
93  1.8874368E7     8.912896E7      0.0
94  5.519704064E9       -1.241513984E9      0.0
95  1.92937984E8        5.0331648E8     0.0
96  1.879048192E10      -2.6306674688E10    0.0
97  4.69762048E9        4.966055936E9       0.0
98  4.7244640256E10     -4.4023414784E10    0.0
99  2.147483648E9       2.147483648E9       0.0

Original comment by os...@united-forum.de on 2 Oct 2013 at 10:48

GoogleCodeExporter commented 9 years ago 
Code snippet to reproduce (both are now based on the same random values, 
hopefully since Colt also uses row-major-ordering):

for(int d = 1; d < 100; d++) {
  // Dimension
  System.out.print(d + "\t");

  double[][] values = new double[d][d];
  Random RNG = new Random();

  for (int i = 0; i < values.length; i++) {
    for (int j = 0; j < values[i].length; j++) {
      values[i][j] = RNG.nextDouble();
    }
  }

  // Error using Colt
  DenseDoubleMatrix2D matrixCern = new DenseDoubleMatrix2D(values);
  Algebra algebra = new Algebra();
  System.out.print(algebra.det(matrixCern) - algebra.det(algebra.transpose(matrixCern)) + "\t");

  // Error using EJML ...
  Random random = new Random();
  DenseMatrix64F matrix = new DenseMatrix64F(values);
  System.out.print(CommonOps.det(matrix) - CommonOps.det(CommonOps.transpose(matrix, null)) + "\t");
  // ... for det(A) - det((A^t)^t), to differentiate whether det() or ^t seems to be erroneous.
  System.out.println(CommonOps.det(matrix) - CommonOps.det(CommonOps.transpose(CommonOps.transpose(matrix, null), null)));
}

Original comment by os...@united-forum.de on 2 Oct 2013 at 10:49

GoogleCodeExporter commented 9 years ago 
It certainly did look weird.  The problem might just be because the determinant 
is very large.  Here is what the error looks like when it is rescaled by the 
magnitude:

det(a)-det(a')
ans = -520

(det(a)-det(a'))/det(a)
ans =  8.8078e-15

Original comment by peter.ab...@gmail.com on 2 Oct 2013 at 11:09

GoogleCodeExporter commented 9 years ago 
Below are the results for the relative error (rescaled by det(A)), which seems 
to be "stable".  I am still quite surprised by the fact that the absolute error 
of the algorithm increases so quickly with regards to the dimension.

Thanks again for the fast reply :-)

Results for (det(A) - det(A^t)) / det(A)
1   0.0         0.0
2   -0.0            -0.0
3   -0.0            -0.0
4   1.824441488823891E-16   1.824441488823891E-16
5   1.206920119610065E-15   -4.023067065366889E-16
6   0.0         -1.3428689459092415E-16
7   -2.8129604068761855E-15 2.163815697597064E-16
8   -8.604406468453307E-16  -2.3662117788246626E-15
9   -2.7105802646278556E-15 -5.70648476763758E-16
10  -0.0            -1.1863194498221541E-16
11  -6.074482570322987E-16  -8.099310093763981E-16
12  -1.729802377919336E-16  -0.0
13  -4.627622848423366E-15  3.3054448917309567E-15
14  0.0         -1.4352639892649032E-15
15  1.3567507497480756E-13  7.672226889947382E-14
16  3.2231009313347807E-16  -1.9338605588008712E-15
17  -3.8228562721225634E-15 -1.7521424580561732E-15
18  9.997976599385661E-16   1.888506690995067E-15
19  -6.247301478393343E-16  -5.497625300986159E-15
20  -3.7817716192992067E-14 1.3789885462637529E-14
21  1.3257632294239314E-15  1.4914836331019226E-15
22  -3.391086168759087E-16  2.2042060096934056E-15
23  -1.1932444388975158E-15 -4.772977755590059E-16
24  -1.3527729664966746E-14 2.0419214588628925E-15
25  5.670520777691471E-13   2.3057251956158276E-13
26  -9.728603632988551E-16  3.729298059312255E-15
27  7.601173118479366E-15   -7.812316816214996E-15
28  3.953139494077823E-15   2.888832707210715E-15
29  -5.6842554005956695E-15 -4.440824531715365E-15
30  -3.605704636428275E-15  -3.605704636428263E-15
31  -4.3905315117264786E-15 -9.220116174625617E-15
32  8.274404030208207E-15   1.155294147613977E-14
33  -3.526813576366259E-15  -8.942991568643024E-15
34  -2.3403138353939926E-16 -1.1701569176969995E-16
35  -4.109194601518285E-15  -2.301148976850253E-15
36  -1.2902167445123622E-15 1.4075091758316651E-15
37  9.167595566845657E-16   -9.167595566845628E-16
38  -6.021428015681001E-15  1.254464169933536E-15
39  1.4666345060552747E-15  -1.0083112229130145E-14
40  5.637349574566548E-15   -4.3998825947836655E-15
41  6.706256346968483E-16   -1.9112830588860036E-14
42  2.6484328156639312E-14  1.6420283457116783E-14
43  -2.369582275869506E-15  -2.1541657053359203E-16
44  4.370918348461875E-15   -6.712481749423625E-15
45  -1.3610540851897579E-15 -8.71074614521448E-15
46  1.1221363621600939E-14  4.98727272071156E-15
47  6.538054155955505E-15   2.124867600685546E-15
48  -1.4405870482795975E-15 -3.781541001733941E-15
49  2.9575547105016463E-14  2.160700591609985E-14
50  5.371192946031127E-16   6.08735200550192E-15
51  -6.403539546510777E-15  -2.022170383108658E-15
52  2.5915609515594963E-16  7.385948711944562E-15
53  -1.734728425344786E-14  -7.059941265938121E-15
54  -1.9534733138010694E-14 1.7063471716936757E-14
55  3.7159207638150917E-14  -5.955209259040166E-14
56  -1.6499754968925105E-14 -3.4020113337989656E-15
57  -1.1276188001243333E-14 -3.3476183128690753E-15
58  2.2142936535296792E-16  -9.078603979471747E-15
59  2.6195912977776717E-15  -1.007535114529887E-14
60  -7.218352908713267E-15  8.478700241980608E-15
61  4.862190273069208E-14   1.435776186517924E-13
62  1.2988855662485257E-14  -8.962310407114962E-15
63  -5.658997373946862E-15  -1.5433629201673266E-15
64  2.431250334383285E-13   -4.76682078678957E-13
65  7.074369030257145E-15   -5.264646720191408E-15
66  -2.4194913067697545E-14 -1.9606222658306205E-14
67  2.947617799949015E-15   -1.4738088999745222E-14
68  2.4286851190253887E-14  1.6874630372448863E-14
69  -4.445551216908515E-15  3.810472471635848E-15
70  3.3609724209992E-15 -2.868029799252456E-14
71  -1.3589070386757636E-15 9.706478847684006E-16
72  -7.166334982735422E-15  8.745357945032958E-15
73  8.492006754479673E-15   -3.569104288114665E-15
74  -2.962006382824995E-15  -3.7594696397394136E-15
75  1.2562495539536512E-14  1.5330503031298824E-14
76  -2.552930253167302E-14  2.3497119245568902E-14
77  3.5304772438571084E-15  -1.922148721655588E-14
78  7.055066079538765E-15   -4.385207739635053E-14
79  8.523075328026016E-15   2.0889890509868004E-14
80  -4.81174680867302E-15   4.277108274375968E-15
81  1.4987386309760915E-14  4.995795436586945E-15
82  -5.788268509573278E-14  4.488552865125873E-14
83  1.8521702981816478E-14  1.438876099413821E-14
84  -1.4199364969082557E-14 -5.409281892983787E-16
85  -0.0            -8.994127563298831E-15
86  -8.059996059903619E-16  3.089665156296395E-15
87  -7.778038821309324E-16  -2.644533199245204E-15
88  -3.744325346365828E-15  -1.5841376465393997E-15
89  1.0332265573379568E-13  -5.942740983380247E-14
90  0.0         1.5999003321621432E-14
91  -4.611493620212044E-14  -1.7877891933431448E-14
92  2.5119513874937508E-14  -3.6075897586347066E-15
93  -2.7657005469684424E-14 -2.228364440700316E-14
94  -2.2271996746778896E-14 1.3622677621816062E-14
95  1.3426279255787078E-14  -2.2057458777364554E-14
96  -1.8108939130105014E-14 -4.527234782526444E-15
97  -2.828897622020963E-14  -4.2365133469878996E-15
98  -2.463523591844831E-15  -1.4429209609376787E-14
99  1.6092441856954997E-16  6.919749998490604E-15

Original comment by os...@united-forum.de on 2 Oct 2013 at 11:41

GoogleCodeExporter commented 9 years ago

Original comment by peter.ab...@gmail.com on 2 Oct 2013 at 11:46