Re: Possible bug in det (was: Possible loss of accuracy)

[louis scott]

Perfect Jordi, that is it. I will submit a report.  

I think the easiest answer is that once it has been LU factored, all we want is det(U) not det(A).  And then the following is correct:

http://hg.savannah.gnu.org/hgweb/octave/file/7f6f0b3f7369/liboctave/array/dMatrix.cc#l1263

recall that for a 3 by 3 case:

C = [ a b c; d e f; g h i];

det(C) = aei + bfg + cdh - ceg - bdi - afh;

but if we now have U = [a' b' c'; 0 e' f'; 0 0 i'];  where prime denotes the transformed coordinates.

then all the terms in det(U)  except the first contain a zero and cancel out.

hence:

if (typ == MatrixType::Lower || typ == MatrixType::Upper)

1264 {

1265 for (octave_idx_type i = 0; i < nc; i++)

1266 retval *= elem (i,i);

1267 }

On Thu, May 16, 2013 at 2:47 PM, Jordi Gutiérrez Hermoso <address@hidden> wrote:

On 16 May 2013 06:13, Marco Caliari <address@hidden> wrote:
> Here I suspect a bug: I get
>
> N = 2100;
> A = rand(N);
> det(A*inv(A))
> ans = Inf
> prod(diag(chol(A*inv(A))))
> ans =  1.00000 % <- correct
> prod(diag([~,U]=lu(A*inv(A))))
> ans =  1.00000 % <- correct
>
> Shouldn't Octave use Cholesky factorization of A*inv(A) (or, if not
> recognized as hermitian, LU factorization) in order to compute its
> determinant?

It is using an LU factorisation. This determinant is computed with
dgetrf, which gives a beautiful LU factorisation. But then... Octave
is doing something funny to actually get the determinant out of that
factorisation. Look at these lines:

     http://hg.savannah.gnu.org/hgweb/octave/file/7f6f0b3f7369/liboctave/array/dMatrix.cc#l1357

Here atmp is the matrix that has been LU-factorised by dgetrf. So it
should just be a simple matter of getting the product of its diagonal,
in retval. But here there is some trickery going on, because retval
here is of det type, and that *= operator is actually overloaded here:

    http://hg.savannah.gnu.org/hgweb/octave/file/7f6f0b3f7369/liboctave/numeric/DET.h#l70

What is that doing? It seems to be doing multiplication by taking
logarithms! The definition of xlog2 is this one:

    http://hg.savannah.gnu.org/hgweb/octave/file/7f6f0b3f7369/liboctave/numeric/lo-mappers.cc#l135

The actual determinant (base_det::value) is then computed as the
result of ldexp () after taking the significand and exponent from
frexp (). But the problem is that here the significand gets very small
and the exponent gets very big, so ldexp () ends up reporting
infinity.

Stepping through the code, this does indeed look like a bug. I don't
understand why it's doing this at all instead of just ordinary
floating point multiplication. I'm guessing it's trying to work around
problems when losing precision due to multiplying numbers that are of
very different magnitude. Please open a bug report, and perhaps
someone else can figure out what the problem is.

- Jordi G. H.