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

 
 

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

Submitter:  Marco Caliari <caliari>
Submitted:  Fri 17 May 2013 01:37:41 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Fri 11 Aug 2017 05:13:39 PM UTC, comment #9: 

Sorry this was ignored for so long.  I pushed the changeset here:

http://hg.savannah.gnu.org/hgweb/octave/rev/6d73cb87591a

John W. Eaton <jwe>
Group administrator
Thu 23 Mar 2017 07:34:13 PM UTC, comment #8: 

I re-encountered this bug. The patch proposed by Andrew solves the problem.

Marco Caliari <caliari>
Group Member
Fri 01 Jul 2016 12:37:32 PM UTC, comment #7: 

The patch looks good to me.

I freshened it to apply to the current default, and rewrote the commit message.

I haven't yet tested it because my system is refusing to build at the moment, but it is so simple I can't imagine it not working.

(file #37650)

Lachlan Andrew <lachlan>
Thu 04 Jun 2015 02:00:08 AM UTC, comment #6: 

This bug had a patch attached bug doesn't seem to have been reviewed or commented on. Updating bug status and release info.

Mike Miller <mtmiller>
Group Member
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>
Group Member
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


A = rand(2100);
det(A*inv(A))


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


prod([1e-200*ones(1,20),1e200*ones(1,20)])
ans =  1.00000


but


prod([1e-200*ones(1,30),1e200*ones(1,30)])
ans = 0


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


function ret = det(A)
if ishermitian(A)
   ret = prod(diag(chol(A)));
else
   ret = prod(diag([~,U]=lu(A)));
end


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

Cheers,

Marco

Marco Caliari <caliari>
Group Member

 

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

Attach Files:
   
   
Comment:
   

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

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by lachlan (Updated the item)
  • -email is unavailable- added by mtmiller (Posted a comment)
  • -email is unavailable- added by jordigh (Posted a comment)
  • -email is unavailable- added by drizzd (Posted a comment)
  • -email is unavailable- added by caliari (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 9 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2017-08-11 jwe StatusPatch Reviewed Fixed
        Open/ClosedOpen Closed
    2017-03-23 caliari StatusPatch Submitted Patch Reviewed
    2016-07-01 lachlan Attached File- Added bug_39014_det_underflow.cset, #37650
    2015-06-04 mtmiller Release3.6.2 dev
        Operating SystemGNU/Linux Any
    2015-06-04 mtmiller Item GroupNone Incorrect Result
        StatusNone Patch Submitted
    2013-05-20 drizzd Attached File- Added det.patch, #28118

    Back to the top

    Powered by Savane 3.13-4448.
    Corresponding source code