Bug: det() gives false zero

[Jaroslav Hajek]

On Thu, Jul 24, 2008 at 4:15 PM, Charlie Dyson <charlie at charliedyson.net> wrote:
> I've been using octave to try to find zeros of the determinants of
> various matrices depending upon some parameter, but have found that
> when the entries in such matrices have wildly different magnitudes
> det() may claim that the matrix is singular when it is not. Using
> Octave 3.0.1 running on XP.
>
> Unfortunately I'll be going on holiday in a few hours so probably
> won't be able to get back to you for a week or so.
>
> Here's some code that reproduces the problem, along with a workaround
> I've been using to cope with it.
>
> Incidentally, I couldn't see any email address at
> http://www.gnu.org/software/octave/bugs.html, though I was able to
> find one in the pdf manual.
>
> Many thanks,
>
> Charlie Dyson
>
> -- CUT --
> D = [ -1.001254707117086 0.00125470711708553 1 0 0 0
>  0.00125470711708553 -1.001254707117086 0 1 0 0
>  0 0 -1.050021149616144 -280214911531872 1 0
>  0 0 -5.590061183230268e-021 -0.9499793732213449 0 1
>  1 0 0 0 0 0
>  0 0 0 0 0 1 ];
>
> disp("Compute determinant with det function:")
> disp( det(D) );
> disp("Should not have been zero!");
> disp( "" );
>
> disp("Now try LU decomposition");
> [l,u,p] = lu( D );
> disp("Should all be non-zero:");
> disp( det(l) );
> disp( det(u) );
> disp( det(p) );
> disp("Again we have a false zero");
> disp("");
>
> disp("Now compute det of L, U by product of diagonal");
> detL = 1;
> for n = 1:rows(D)
>        detL = detL*l(n,n);
> endfor
> detU = 1;
> for n = 1:rows(D)
>        detU = detU*u(n,n);
> endfor
> disp("Det of l,u,p given by:");
> disp( detL );
> disp( detU );
> disp( det(p) );
>
> disp("Correct value of determinant is:");
> disp( detL*detU*det(p) )
> _______________________________________________
> Bug-octave mailing list
> Bug-octave at octave.org
> https://www.cae.wisc.edu/mailman/listinfo/bug-octave
>

Looking at src/DLD-FUNCTIONS/det.cc, it seems that this is a feature -
if a numerically singular matrix is detected, zero is returned.
Perhaps John can comment on why he made this choice.

-- 
RNDr. Jaroslav Hajek
computing expert
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz