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 May 17 13:37:41 2013  
Category: Octave FunctionSeverity: 3 - Normal
Priority: 5 - NormalItem Group: Incorrect Result
Status: Patch ReviewedAssigned to: None
Originator Name: Open/Closed: Open
Release: devOperating System: Any

Add a New Comment (Rich MarkupRich Markup):

You are not logged in

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


(Jump to the original submission Jump to the original submission)

Thu Mar 23 19:34:13 2017, comment #8:

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

Marco Caliari <caliari>
Project Member
Fri Jul 1 12:37:32 2016, 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>
Project Member
Thu Jun 4 02:00:08 2015, 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>
Project Administrator
Mon May 20 21:07:18 2013, 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 May 20 13:48:52 2013, 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 May 18 11:57:08 2013, 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 May 18 11:03:09 2013, 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 May 18 10:49:51 2013, 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}
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 May 17 13:37:41 2013, 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


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).



Marco Caliari <caliari>
Project Member


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

Attach File(s):

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 lachlan (Updated the item)
  • -unavailable- added by mtmiller (Posted a comment)
  • -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):



    Follow 7 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Thu Mar 23 19:34:13 2017caliariStatusPatch Submitted=>Patch Reviewed
    Fri Jul 1 12:37:32 2016lachlanAttached File-=>Added bug_39014_det_underflow.cset, #37650
    Thu Jun 4 02:00:08 2015mtmillerRelease3.6.2=>dev
      Operating SystemGNU/Linux=>Any
    Thu Jun 4 02:00:07 2015mtmillerItem GroupNone=>Incorrect Result
      StatusNone=>Patch Submitted
    Mon May 20 13:48:52 2013drizzdAttached File-=>Added det.patch, #28118

    Back to the top

    Powered by Savane 3.1-cleanup1