Forum: Ruby Is this a floating point precision problem?

Announcement (2017-05-07): www.ruby-forum.com is now read-only since I unfortunately do not have the time to support and maintain the forum any more. Please see rubyonrails.org/community and ruby-lang.org/en/community for other Rails- und Ruby-related community platforms.
Todd S. (Guest)
on 2006-01-30 22:09
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?
Eero S. (Guest)
on 2006-01-30 23:29
(Received via mailing list)
On 2006.01.31 05:09, Todd S. wrote:
> will invert to this...
> 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 C. (Guest)
on 2006-01-31 00:06
(Received via mailing list)
On Jan 30, 2006, at 4:28 PM, Eero S. wrote:

>>     -0.998844486916645      0.0480592442327619 0.0
>> multiplyng by the inverse matrix and then back again by the matrix
>
>> 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 T. (Guest)
on 2006-01-31 01:06
(Received via mailing list)
> 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. (Guest)
on 2006-01-31 02:42
>
> 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]]
Todd S. (Guest)
on 2006-01-31 20:28
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?
M. Edward (Ed) Borasky (Guest)
on 2006-02-01 17:32
(Received via mailing list)
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:
> will invert to this...
> 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://linuxcapacityplanning.com
Jacob F. (Guest)
on 2006-02-01 17:59
(Received via mailing list)
On 2/1/06, M. Edward (Ed) Borasky <removed_email_address@domain.invalid> wrote:
>
> 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 F.

[1]
http://www.ruby-doc.org/stdlib/libdoc/matrix/rdoc/...
Zinsser (Guest)
on 2006-02-02 00:40
(Received via mailing list)
Todd S. wrote:

> -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 ...
Zinsser (Guest)
on 2006-02-02 01:20
(Received via mailing list)
> 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 (Guest)
on 2006-02-03 08:03
(Received via mailing list)
Yeah ... I think if you use Mathn you'll have better results ...
probably slower though

Zinsser wrote:
> [-0.974313737622412, -0.225194894880513, -2.52327742604428e-17]]
>     for k in 0..size
>
>     for k in 0..size
>
>
>
>

--
M. Edward (Ed) Borasky

http://linuxcapacityplanning.com
This topic is locked and can not be replied to.