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

[louis scott]

Sergei, 

Thank you.  The rand function calls a Mersenne Twister, so the odds of getting a column to be a linear combination of other columns is tiny.  But it is true that the condition number rises with size, and that the issue resides in inv(A), but the numbers below would suggest that it is reasonable to invert:

N=[10 100 500 1000 2000];

for ith=1:10

     for nth=1:5

        tt(nth,ith) = cond(rand(N(nth)));

     end

end

mean(tt')

1.6747e+02   1.5101e+04   3.6118e+04   3.1278e+05   5.6422e+05

Best regards, Louis

Ah Marco, well done. It is the off diagonals that are noisy.

A = rand(2200);

det(A*inv(A))

        ans = Inf

prod(diag(A*inv(A)))

        ans =  1.00000

So consider this:

B= A*inv(A); % this should be the identity matrix

[~,U] = lu(B);

for ith=1:2200  %focus on the off diagonals

   B(ith,ith)=0;

   U(ith,ith)=0;

end

plot(max(B))  % but the off diagonals are not zero, rather it is noise on the order of 1e-12

plot(max(B-B')) % alternatively asking if it is Hermitian and no - the noise again.

plot(max(U)) % noisy still.

Now for a simple 3 by 3

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

det(C) = aei + bfg + cdh - ceg - bdi - afh;  the last terms contain a 1 in the product from the diagonal

Combining terms:

det(C) = 1 + ((N-1)*O(1e-12^N)) - (N*O(1e-12^(N-1))) 

On Thu, May 16, 2013 at 11:13 AM, Marco Caliari <address@hidden> wrote:

Secondly, this makes use of ATLAS' multi core capabilities:
% getting to 1 the hard way
N = 2100; A = rand(N);
tic, det(A*inv(A)), toc
ans = Inf
Elapsed time is 1.6111 seconds.

N = 2000; A = rand(N);
tic, det(A*inv(A)), toc
ans =  1.0000
Elapsed time is 1.436 seconds.

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?

Marco