Path: news.mathworks.com!not-for-mail
From: "Bruno Luong" <b.luong@fogale.findmycountry>
Newsgroups: comp.soft-sys.matlab
Subject: Re: An interesting numerical problem
Date: Thu, 22 Oct 2009 17:59:07 +0000 (UTC)
Organization: FOGALE nanotech
Lines: 54
Message-ID: <hbq6hb$5me$1@fred.mathworks.com>
References: <hbq4lq$4m1$1@fred.mathworks.com>
Reply-To: "Bruno Luong" <b.luong@fogale.findmycountry>
NNTP-Posting-Host: webapp-05-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1256234347 5838 172.30.248.35 (22 Oct 2009 17:59:07 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Thu, 22 Oct 2009 17:59:07 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 390839
Xref: news.mathworks.com comp.soft-sys.matlab:579400


"Yi Cao" <y.cao@cranfield.ac.uk> wrote in message <hbq4lq$4m1$1@fred.mathworks.com>...
> I have met an interesting numerical problem:
> 
> Let's define a function first:
> 
> function P=polyconvert(n)
> P = eye(n);
> for k=1:n
>      P(k,1:k) = fliplr(poly(ones(1,k-1)));
> end
> 
> Then use the function as follows:
> 
> P19 = polyconvert(19);
> 
> This is a low triangular matrix with all diagonal elements equal to 1. Therefore, the determinant of the matrix should be 1, as well any product of P19 and P19'. However, in Matlab 2009b, I got the following results:
> 
> det(P19)        % I got 1, it is fine.
> det(P19'*P19) % I also got 1.
> det(P19*P19') % I got -8!
> 
> It is very strange, how can det(P19*P19') equal to -8? All elements of P19 are integers, hence should not have any error here. The maximum elements of P19*P19' is only about 1e10. Although it is slightly over the 2^32, but that should not be a problem for double floating point calculation, isn't it?


>Anyone can give a clue what is going wrong here?
> 

cond(P19*P19') is larger than 10^20. Eigen-space are highly parallel. The fact that P19 contains integer does not matter IMO because DET calculation using LU factorization.

>> det(P19*P19')

ans =

    -8

>> det(fliplr(P19*P19'))

ans =

    -1

>> det(flipud(P19*P19'))

ans =

     4

>> det(fliplr(flipud(P19*P19')))

ans =

     1

>>