[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Octave-bug-tracker] [bug #39014] Wrong determinant for some (large) mat
From: |
Clemens Buchacher |
Subject: |
[Octave-bug-tracker] [bug #39014] Wrong determinant for some (large) matrices |
Date: |
Sat, 18 May 2013 10:49:51 +0000 |
User-agent: |
Mozilla/5.0 (X11; Linux x86_64) AppleWebKit/537.31 (KHTML, like Gecko) Chrome/26.0.1410.63 Safari/537.31 |
Follow-up Comment #1, bug #39014 (project octave):
For the determinant we make some effort to deal with large intermediate
results. We normalize factors using frexp and keep track of the exponent
separately, such that we are not limited by the double precision exponent (+/-
1023). Let's assume all factors are close to 1. For those smaller than 1,
frexp returns an exponent of 0 and leaves the mantissa untouched. For those
equal to or larger than 1, frexp returns an exponent of 1 and divides the
mantissa by 2. So for those latter factors we have a mantissa close to 0.5,
and we keep multiplying with 0.5 until we do hit the double precision limit in
the mantissa:
> det(matrix_type(full(diag(ones(1, 2000)-1e-10)), "full"))
ans = 1.00000
> det(matrix_type(full(diag(ones(1, 2000))), "full"))
ans = NaN
Breakpoint 6, Matrix::determinant (address@hidden,
mattype=..., address@hidden: 0, address@hidden: 1,
address@hidden) at ../../liboctave/array/dMatrix.cc:1362
1362 retval *= (ipvt(i) != (i+1)) ? -c : c;
3: i = 1066
2: retval = {c2 = 6.3240402667679558e-322, e2 = 1067}
(gdb)
Will ignore next 9 crossings of breakpoint 6. Continuing.
Breakpoint 6, Matrix::determinant (address@hidden,
mattype=..., address@hidden: 0, address@hidden: 1,
address@hidden) at ../../liboctave/array/dMatrix.cc:1362
1362 retval *= (ipvt(i) != (i+1)) ? -c : c;
3: i = 1076
2: retval = {c2 = 0, e2 = 1077}
I wonder why the limit on my machine appears to be 2^-1074 ~ 1e-325.
So considering the effort we make already, it seems that we should be doing
better here. How about renormalizing retval every 100 iterations or so?
_______________________________________________________
Reply to this item at:
<http://savannah.gnu.org/bugs/?39014>
_______________________________________________
Message sent via/by Savannah
http://savannah.gnu.org/