bugGNU Octave - Bugs: bug #39014, Wrong determinant for some (large)...

 
 

bug #39014: Wrong determinant for some (large) matrices

Submitted by:  Marco Caliari <caliari>
Submitted on:  Fri 17 May 2013 01:37:41 PM UTC  
 
Category: Octave FunctionSeverity: 3 - Normal
Priority: 5 - NormalItem Group: None
Status: NoneAssigned to: None
Originator Name: Open/Closed: Open
Release: 3.6.2Operating System: GNU/Linux

Add a New Comment (Rich MarkupRich Markup):
   

You are not logged in

Please log in, so followups can be emailed to you.

 

Mon 20 May 2013 09:07:18 PM UTC, comment #5:

Clemens Buchacher wrote:

> on yet another octave mailing list.


I am not sure if you're suggesting that we have too many Octave mailing lists, but if so, you should rejoice about the situation getting better. We had three, help, dev, and Octave-Forge, but we recently folded Octave-Forge into the dev list.

There are now only two Octave lists: help and dev.

Jordi GutiƩrrez Hermoso <jordigh>
Project Administrator
Mon 20 May 2013 01:48:52 PM UTC, comment #4:

Attaching det.patch with the following description:

normalize multiplication result instead of operand

The operand of a multiplication with base_det is normalized to c * 2^e,
such that 0.5 <= c < 1. The mantissa is multiplied by c, and the
exponent is updated separately. After each multiplication, the mantissa
is therefore smaller by a up to a factor 2. Despite the normalization of
the operands, the mantissa can therefore become very small. Eventually,
each multiplication suffers from severe loss of precision and the
mantissa can even get rounded to 0.

Instead, normalize the mantissa after multiplication with the operand.
This makes the operation insensitive to the number of operands.

The change may degrade accuracy in case of operands close to the maximum
or minimum representable floating point numbers, but that is not a
problem base_det tries to solve.

Bug report: https://savannah.gnu.org/bugs/?39014
Reported-by: Marco Caliari <marco.caliari@univr.it>
Fix-proposed-by: Stephen Montgomery-Smith <stephen@missouri.edu>

(file #28118)

Clemens Buchacher <drizzd>
Sat 18 May 2013 11:57:08 AM UTC, comment #3:

I see now that this was not only reported, but also analyzed and practically solved already on yet another octave mailing list. Oh well, I guess I should consider this a good exercise.

Clemens Buchacher <drizzd>
Sat 18 May 2013 11:03:09 AM UTC, comment #2:

Err, I veered off course a bit here. But I suspect the original issue is quite similar:

octave-3.7.3+:52> det(matrix_type(full(diag(ones(1,2000)+1e-10)), "full"))
ans = Inf

In this case, we do not end up with c2=0, but we end up with a very small mantissa, which due to loss of precision does not grow as quickly as the exponent. So the final result ends up as retval = c2 * 2^e2 = 2^-1074 * 2^2000, which I suppose gets converted to Inf somewhere (maybe we could do better there too, and get 2^928, but that would not save our day).

Breakpoint 6, Matrix::determinant (this=this@entry=0x7fffffffd2e0,
mattype=..., info=@0x7fffffffd250: 0, rcon=@0x7fffffffd280: 1,
calc_cond=calc_cond@entry=1) at ../../liboctave/array/dMatrix.cc:1362
1362 retval *= (ipvt(i) != (i+1)) ? -c : c;
4: c = 1.0000000001
3: i = 1821
2: retval = {c2 = 4.9406564584124654e-324, e2 = 1822}

Clemens Buchacher <drizzd>
Sat 18 May 2013 10:49:51 AM UTC, comment #1:

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 (this=this@entry=0x7fffffffd2e0,
mattype=..., info=@0x7fffffffd250: 0, rcon=@0x7fffffffd280: 1,
calc_cond=calc_cond@entry=1) 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 (this=this@entry=0x7fffffffd2e0,
mattype=..., info=@0x7fffffffd250: 0, rcon=@0x7fffffffd280: 1,
calc_cond=calc_cond@entry=1) 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?

Clemens Buchacher <drizzd>
Fri 17 May 2013 01:37:41 PM UTC, original submission:

Dear all,

as reported in the thread https://mailman.cae.wisc.edu/pipermail/help-octave/2013-May/058753.html, the following code

quite often gives Inf instead of 1. The problem seems to be in the way the product of the diagonal elements of the R (or U) factor of the matrix is done. I don't know if prod is sufficiently robust with respect to the product of numbers very different in magnitude. I see

but

Suppose prod is fixed if needed. Isn't then it possible to rewrite det function with a simple m file basically doing

By the way, the actual implementation of prod allows the script above to work fine on A*inv(A).

Cheers,

Marco

Marco Caliari <caliari>
Project Member

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach File(s):
   
   
Comment:
   

Attached Files
file #28118:  det.patch added by drizzd (2kB - text/x-patch)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -unavailable- added by jordigh (Posted a comment)
  • -unavailable- added by drizzd (Posted a comment)
  • -unavailable- added by caliari (Submitted the item)
  •  

    Do you think this task is very important?
    If so, you can click here to add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.

     

    Please enter the title of George Orwell's famous dystopian book (it's a date):

     

     

    Follows 1 latest change.

    Date Changed By Updated Field Previous Value => Replaced By
    Mon 20 May 2013 01:48:52 PM UTCdrizzdAttached File-=>Added det.patch, #28118

    Back to the top


    Powered by Savane 3.1-cleanup