Is this a floating point precision problem?

ruby-talk
1

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-forum.com/.

2

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

···

On 2006.01.31 05:09, Todd S. wrote:

3

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]]

4

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://linuxcapacityplanning.com

5

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.

···

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

6

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...

a = Matrix[[-0.0108226964530337, -0.224934679233174,

-0.974313737622412], [0.0468247818757306, 0.973187905351297,
-0.225194894880513], [-0.998844486916645, 0.0480592442327619, 0.0]]

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-forum.com/\.

7

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.html#M001009

···

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

8

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-forum.com/\.

9

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...

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

> 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 ...

10

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

11

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