MATLAB Central Newsreader - An interesting numerical problem

An interesting numerical problem

Yi Cao — Thu, 22 Oct 2009 17:27:22 -0400

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?

Thanks for help.

Yi

Re: An interesting numerical problem

Matt Fetterman — Thu, 22 Oct 2009 17:48:19 -0400

"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?
>
> Thanks for help.
>
> Yi

What happens with N less than 19? If it suddenly doesn't work at 19, that does suggest the memory size? Regards! Matt

Re: An interesting numerical problem

Bruno Luong — Thu, 22 Oct 2009 17:59:07 -0400

"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

>>

Re: An interesting numerical problem

Yi Cao — Thu, 22 Oct 2009 19:18:21 -0400

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hbq6hb$5me$1@fred.mathworks.com>...
>
> 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.
>

For L*U = P*P', theoretically, L should be P and U should be P' because P is low triangular and P' is upper triangular. Clearly, LU factorization in Matlab is not so smart.

Yi

Re: An interesting numerical problem

someone — Thu, 22 Oct 2009 19:35:23 -0400

"Yi Cao" <y.cao@cranfield.ac.uk> wrote in message <hbqb5s$9st$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hbq6hb$5me$1@fred.mathworks.com>...
> >
> > 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.
> >
>
> For L*U = P*P', theoretically, L should be P and U should be P' because P is low triangular and P' is upper triangular. Clearly, LU factorization in Matlab is not so smart.
>
> Yi

IMO, its not so much that "Matlab is not so smart" as it is that due to
floating point arithmetic, numerical problems and theoretical problems
don't always agree. In fact, you can even find "simple" cases where

(a + b) + c ~= a + (b + c)