From: "M. Edward (Ed) Borasky" <znmeb@...>
Date: 2006-02-03T15:02:55+09:00
Subject: Re: Is this a floating point precision problem?

Yeah ... I think if you use Mathn you'll have better results ... 
probably slower though

Zinsser wrote:
>> Could this be architecture related?  I've tried this both on a pc and 
>> a mac and get the same result.  What are you running on?
>>     
>
> Hi Todd,
>
> I looked into matrix.rb and found the reason for your problem.
> BTW, I also get the wrong result when calling a.inv:
>
> Matrix[[-8.0, 0.0, -0.998844486916645], [0.0, 1.0625, 0.0480592442327619],
> [-0.974313737622412, -0.225194894880513, -2.52327742604428e-17]]
>
> It is a floating point precision problem in matrix.rb, or, 
> to be more precise, the problem is that the algorithm 
> for inverting the matrix is not very sophisticated.
>
>
> def inverse_from(src)
>     size = row_size - 1
>     a = src.to_a
>     
>     for k in 0..size
>       if (akk = a[k][k]) == 0
>         i = k
>         begin
>           Matrix.Raise ErrNotRegular if (i += 1) > size
>         end while a[i][k] == 0
>         a[i], a[k] = a[k], a[i]
>         @rows[i], @rows[k] = @rows[k], @rows[i]
>         akk = a[k][k]
>       end
> ....
>
> For your matrix, a[1][1] evaluates to -2.33146835171283e-15,
> which is close to zero, but does not get caught by 
> "if (akk = a[k][k]) == 0". If you replace this version with 
>
> def inverse_from(src)
>     size = row_size - 1
>     a = src.to_a
>
>     puts a
>     
>     for k in 0..size
>       if (akk = a[k][k]) > -1.0e-12 && akk < 1.0e-12
>         i = k
>         begin
>           Matrix.Raise ErrNotRegular if (i += 1) > size
>         end while a[i][k] == 0
>         a[i], a[k] = a[k], a[i]
>         @rows[i], @rows[k] = @rows[k], @rows[i]
>         akk = a[k][k]
>       end
> ....
>
> you get the correct result. Perhaps someone could implement 
> this algorithm with a better pivoting strategy?
>
> Your matrix seems to be orthonormal, aka a 3-D rotation matrix.
> For this kind of matrix, the inverse is identical to the transpose,
> which is also much faster to compute.
>
> Timo Zinsser
>
>
>
>
>   

-- 
M. Edward (Ed) Borasky

http://linuxcapacityplanning.com