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?

on 2006-01-30 22:09

on 2006-01-30 23:29

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

on 2006-01-31 00:06

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.

on 2006-01-31 01:06

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

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

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?

on 2006-02-01 17:32

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

on 2006-02-01 17:59

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

on 2006-02-02 00:40

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

on 2006-02-02 01:20

> 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

on 2006-02-03 08:03

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