Is this a floating point precision problem?

[Todd S.]

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

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]

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.


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]

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

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

 

[Zinsser]

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

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?

 

[Zinsser]

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[k] == 0
a, a[k] = a[k], a
@rows, @rows[k] = @rows[k], @rows
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[k] == 0
a, a[k] = a[k], a
@rows, @rows[k] = @rows[k], @rows
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]

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.

 

[Jacob Fugal]

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.html#M=
001009

 

[M. Edward (Ed) Borasky]

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

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[k] == 0
a, a[k] = a[k], a
@rows, @rows[k] = @rows[k], @rows
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[k] == 0
a, a[k] = a[k], a
@rows, @rows[k] = @rows[k], @rows
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