[lnkForumImage]
TotalShareware - Download Free Software
Confronta i prezzi di migliaia di prodotti.
Asp Forum
 Home | Login | Register | Search 


 

Forums >

comp.lang.ruby

Is this a floating point precision problem?

Todd S.

1/30/2006 8:10:00 PM

I'm using the Matrix module (matrix.rb) to help place a vertex into
local coordinates and then back into world coordinates. Under certain
situations the transformation to and from local space corrupts the
data...

This matrix...
-0.0108226964530337 -0.224934679233174 -0.974313737622412
0.0468247818757306 0.973187905351297 -0.225194894880513
-0.998844486916645 0.0480592442327619 0.0

will invert to this...
64.0 0.0
-0.998844486916645
0.0 1.0
0.0480592442327619
-0.974313737622412 -0.225194894880513 -2.40312135813741e-18


starting with vertex: 0.015525, -0.28474399, 0.53221899
multiplyng by the inverse matrix and then back again by the matrix
itself will not yield the same number as it should but instead
returns... 0.00555723611485535, -0.24161810434515, -0.4739174731596

So am I seeing strange behavior due to the extreamly small (e-18) number
in the matrix?

If so, how can I set the precision so that I don't see this error?

--
Posted via http://www.ruby-....
10 Answers

ES

1/30/2006 9:28:00 PM

0

On 2006.01.31 05:09, Todd S. wrote:
> I'm using the Matrix module (matrix.rb) to help place a vertex into
> local coordinates and then back into world coordinates. Under certain
> situations the transformation to and from local space corrupts the
> data...
>
> This matrix...
> -0.0108226964530337 -0.224934679233174 -0.974313737622412
> 0.0468247818757306 0.973187905351297 -0.225194894880513
> -0.998844486916645 0.0480592442327619 0.0
>
> will invert to this...
> 64.0 0.0
> -0.998844486916645
> 0.0 1.0
> 0.0480592442327619
> -0.974313737622412 -0.225194894880513 -2.40312135813741e-18
>
>
> starting with vertex: 0.015525, -0.28474399, 0.53221899
> multiplyng by the inverse matrix and then back again by the matrix
> itself will not yield the same number as it should but instead
> returns... 0.00555723611485535, -0.24161810434515, -0.4739174731596
>
> So am I seeing strange behavior due to the extreamly small (e-18) number
> in the matrix?

Without pretending to have any idea about mathematical operations,
it is quite likely that floating-pointedness is causing the problem.
You cannot really even count on 1.8 and 2.1 - 0.3 being the same.

> If so, how can I set the precision so that I don't see this error?

BigDecimal in the stdlib is, as far as I know, lossless but of
course somewhat slower if you are doing anything very intense.



E

Logan Capaldo

1/30/2006 10:05:00 PM

0


On Jan 30, 2006, at 4:28 PM, Eero Saynatkari wrote:

> On 2006.01.31 05:09, Todd S. wrote:
>> I'm using the Matrix module (matrix.rb) to help place a vertex into
>> local coordinates and then back into world coordinates. Under
>> certain
>> situations the transformation to and from local space corrupts the
>> data...
>>
>> This matrix...
>> -0.0108226964530337 -0.224934679233174 -0.974313737622412
>> 0.0468247818757306 0.973187905351297 -0.225194894880513
>> -0.998844486916645 0.0480592442327619 0.0
>>
>> will invert to this...
>> 64.0 0.0
>> -0.998844486916645
>> 0.0 1.0
>> 0.0480592442327619
>> -0.974313737622412 -0.225194894880513 -2.40312135813741e-18
>>
>>
>> starting with vertex: 0.015525, -0.28474399, 0.53221899
>> multiplyng by the inverse matrix and then back again by the matrix
>> itself will not yield the same number as it should but instead
>> returns... 0.00555723611485535, -0.24161810434515, -0.4739174731596
>>
>> So am I seeing strange behavior due to the extreamly small (e-18)
>> number
>> in the matrix?
>
> Without pretending to have any idea about mathematical operations,
> it is quite likely that floating-pointedness is causing the problem.
> You cannot really even count on 1.8 and 2.1 - 0.3 being the same.
>
>> If so, how can I set the precision so that I don't see this error?
>
> BigDecimal in the stdlib is, as far as I know, lossless but of
> course somewhat slower if you are doing anything very intense.
>
>
>
> E
>

Alternatively, if you know how precise you need it to be, and its not
too precise, and speed is important you can use fixed-point math
instead of floating point.

Caleb Tennis

1/30/2006 11:05:00 PM

0

>
> This matrix...
> -0.0108226964530337 -0.224934679233174 -0.974313737622412
> 0.0468247818757306 0.973187905351297 -0.225194894880513
> -0.998844486916645 0.0480592442327619 0.0
>
> will invert to this...
> 64.0 0.0
> -0.998844486916645
> 0.0 1.0
> 0.0480592442327619
> -0.974313737622412 -0.225194894880513 -2.40312135813741e-18
>
>

I don't get the same inverse; maybe that's a good place to look?

irb(main):031:0> a.inv
=> Matrix[[-0.0108184814453125, 0.0468254089355469,
-0.998844487934061], [-0.224935054779053, 0.973188042640686,
0.0480592423711387], [-0.974313745084977, -0.22519489645261,
-1.18537417594384e-11]]


With my inverse above, v*a*ainv looks good:

irb(main):033:0> v*a*ainv
=> Matrix[[0.0155250005999241, -5.8380675554276e-10,
-6.04039249192278e-21], [4.76057619505643e-08, -0.284744036326806,
-1.09843285920372e-18], [-2.25029812272581e-06,
-3.29933630031226e-07, 0.53221899]]

Todd S.

1/31/2006 12:43:00 AM

0


>
> I don't get the same inverse; maybe that's a good place to look?
>

That's strange. I ran it again and this time came up with a different
number...

irb> a = Matrix[[-0.0108226964530337, -0.224934679233174,
-0.974313737622412], [0.0468247818757306, 0.973187905351297,
-0.225194894880513], [-0.998844486916645, 0.0480592442327619, 0.0]]

irb> puts a.inverse

Matrix[[8.0, 2.0, -0.998844486916645], [-0.25, 1.0, 0.0480592442327619],
[-0.974313737622412, -0.225194894880513, -2.16280922232366e-17]]

--
Posted via http://www.ruby-....

Zinsser

1/31/2006 10:27:00 AM

0

Todd S. wrote:

>
>>
>> I don't get the same inverse; maybe that's a good place to look?
>>
>
> That's strange. I ran it again and this time came up with a different
> number...
>
> irb> a = Matrix[[-0.0108226964530337, -0.224934679233174,
> -0.974313737622412], [0.0468247818757306, 0.973187905351297,
> -0.225194894880513], [-0.998844486916645, 0.0480592442327619, 0.0]]
>
> irb> puts a.inverse
>
> Matrix[[8.0, 2.0, -0.998844486916645], [-0.25, 1.0, 0.0480592442327619],
> [-0.974313737622412, -0.225194894880513, -2.16280922232366e-17]]
>

I tried it in C++

a:
-0.0108227 -0.224935 -0.974314
0.0468248 0.973188 -0.225195
-0.998844 0.0480592 0

inverse of a:
-0.0108227 0.0468248 -0.998844
-0.224935 0.973188 0.0480592
-0.974314 -0.225195 -7.06138e-16

a * inverse:
1 2.65358e-17 1.1294e-18
-6.53232e-18 1 -3.77082e-18
-8.77526e-19 9.58502e-18 1

Caleb is using a.inv,
whereas you are using a.inverse.

Perhaps that's your problem ...

Todd S.

1/31/2006 6:29:00 PM

0

Todd S. wrote:
>
>>
>> I don't get the same inverse; maybe that's a good place to look?
>>


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?

--
Posted via http://www.ruby-....

Zinsser

2/1/2006 9:48:00 AM

0

> 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

2/1/2006 3:31:00 PM

0

I don't know what Ruby is doing here, but the R language gives a much
different result for the matrix.

Your input matrix:

> m1
[,1] [,2] [,3]
[1,] -0.01082270 -0.22493468 -0.9743137
[2,] 0.04682478 0.97318791 -0.2251949
[3,] -0.99884449 0.04805924 0.0000000

R's inverse, computed using the routine "ginv"

> library(MASS)
> m3<-ginv(m1)
> m3
[,1] [,2] [,3]
[1,] -0.01082270 0.04682478 -9.988445e-01
[2,] -0.22493468 0.97318791 4.805924e-02
[3,] -0.97431374 -0.22519489 6.223320e-17

Multiplication check:
> m1%*%m3
[,1] [,2] [,3]
[1,] 1.000000e+00 -1.509074e-16 -8.675039e-17
[2,] -1.865234e-16 1.000000e-00 -1.023683e-17
[3,] 6.958291e-17 -1.558541e-18 1.000000e+00

So ... either the Ruby matrix inverter is broken or you matrix is
ill-conditioned enough that the Ruby matrix inverter can't handle it. I
don't have a condition number checker handy, but I can do a singular
value decomposition and look at the numerical rank.

> svd(m1)
$d
[1] 1 1 1

$u
[,1] [,2] [,3]
[1,] -0.01082270 -0.22493468 -9.743137e-01
[2,] 0.04682478 0.97318791 -2.251949e-01
[3,] -0.99884449 0.04805924 6.223320e-17

$v
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 1 0
[3,] 0 0 1

Hmmm ... no zero or small singular values ("$d") so ... looks off the
top of my head like a nice rank 3 well conditioned matrix. So ... best
have a look at the Ruby matrix inverter.


Todd S. wrote:
> I'm using the Matrix module (matrix.rb) to help place a vertex into
> local coordinates and then back into world coordinates. Under certain
> situations the transformation to and from local space corrupts the
> data...
>
> This matrix...
> -0.0108226964530337 -0.224934679233174 -0.974313737622412
> 0.0468247818757306 0.973187905351297 -0.225194894880513
> -0.998844486916645 0.0480592442327619 0.0
>
> will invert to this...
> 64.0 0.0
> -0.998844486916645
> 0.0 1.0
> 0.0480592442327619
> -0.974313737622412 -0.225194894880513 -2.40312135813741e-18
>
>
> starting with vertex: 0.015525, -0.28474399, 0.53221899
> multiplyng by the inverse matrix and then back again by the matrix
> itself will not yield the same number as it should but instead
> returns... 0.00555723611485535, -0.24161810434515, -0.4739174731596
>
> So am I seeing strange behavior due to the extreamly small (e-18) number
> in the matrix?
>
> If so, how can I set the precision so that I don't see this error?
>
>

--
M. Edward (Ed) Borasky

http://linuxcapacitypl...

Jacob Fugal

2/1/2006 3:59:00 PM

0

On 2/1/06, M. Edward (Ed) Borasky <znmeb@cesmail.net> wrote:
> I don't know what Ruby is doing here, but the R language gives a much
> different result for the matrix.
>
> Your input matrix:
>
> > m1
> [,1] [,2] [,3]
> [1,] -0.01082270 -0.22493468 -0.9743137
> [2,] 0.04682478 0.97318791 -0.2251949
> [3,] -0.99884449 0.04805924 0.0000000
>
> R's inverse, computed using the routine "ginv"
>
> > library(MASS)
> > m3<-ginv(m1)
> > m3
> [,1] [,2] [,3]
> [1,] -0.01082270 0.04682478 -9.988445e-01
> [2,] -0.22493468 0.97318791 4.805924e-02
> [3,] -0.97431374 -0.22519489 6.223320e-17

Comparing your inverse and Caleb's, they look pretty close. In only
did a cursory check of the first few digits in each number, but for
that they all matched up except the lower-right number. In that slot
Caleb's was of magnitude e-11 and R's is e-17. I blame that
discrepancy on different approaches and/or convergence criteria
(guessing that it's using a numerical approach).

There's big differences however between Caleb's inverse and Todd's,
same as between yours and Todd's. My only guess is that looking at
Todd's irb session he seems to be using Matrix#inverse while Caleb is
using Matrix#inv. The documentation[1], however, states that
Matrix#inv is just an alias for Matrix#inverse. So I'm not sure what
this is.

Todd, are there any other libraries you're loading which may be
changing the definition of Matrix#inverse?

Jacob Fugal

[1] http://www.ruby-doc.org/stdlib/libdoc/matrix/rdoc/classes/Matrix.ht...

M. Edward (Ed) Borasky

2/3/2006 6:03:00 AM

0

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://linuxcapacitypl...