From: noreply@... Date: 2007-03-19T09:16:02+09:00 Subject: [ ruby-Bugs-9360 ] Matrix inverse problem Bugs item #9360, was opened at 2007-03-18 14:41 You can respond by visiting: http://rubyforge.org/tracker/?func=detail&atid=1698&aid=9360&group_id=426 Category: Standard Library Group: 1.8.6 Status: Open Resolution: None Priority: 3 Submitted By: Stefan M��hlebach (muehlebach) Assigned to: Nobody (None) Summary: Matrix inverse problem Initial Comment: There is a problem with the Matrix.inverse method. Given the following matrix m = Matrix[[0.417787968829298, 0.542037605627772, 0.729142268138943], [0.0, 6.12303176911189e-17, 1.0], [-0.542037605627772, 0.417787968829298, -2.55812900589452e-17]] the result of m.inverse * m is [[1.0, 0.0, 0.650423523448542], [5.55111512312578e-17, 1.0, 0.84385869292008], [0.0, 9.24446373305873e-33, 1.0]] but it should be [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]] or close to this! ---------------------------------------------------------------------- Comment By: Ed Borasky (znmeb) Date: 2007-03-18 16:15 Message: I used R, not Ruby, to do the computations. However, there is an R - Ruby interface called "RSRuby" on RubyForge, so if you're happier in Ruby, you can create your data in Ruby and pass the matrix off to R for computation. For just plain linear equation solving, doing the inverse is usually (almost always) a bad idea. There are lots of alternatives in Ruby, including the Ruby - GNU Scientific Library (GSL) interface, and the built-in equation solver in the "BigDecimal" package. If you have your heart set on pure Ruby, check out the BigDecimal version first. This could actually be a bug in the inverse routine in "Matrix". I looked at it and it looks like it was copied out of an ancient FORTRAN scientific library, and there aren't any comments. It might just be a typo somewhere in the code. Since I don't in general do matrix inverses, and when I do, I usually use someone else's code. If I can find the original FORTRAN (I think I have it around the house somewhere, believe it or not) :) I'll check it. ---------------------------------------------------------------------- Comment By: Stefan M��hlebach (muehlebach) Date: 2007-03-18 15:57 Message: 1. Yes, I did load 'mathn'. 2. I'm trying to solve an equation system. I did exactly the same using Java and the J3D-Library where the results where more accurate. 3. Your results look much better! Where do I get the MASS library? ---------------------------------------------------------------------- Comment By: Ed Borasky (znmeb) Date: 2007-03-18 15:27 Message: Questions: 1. Have you loaded "mathn"? If not, load it and try again. The "Matrix" routines are more accurate with "mathn" loaded. 2. What is the problem you're actually trying to solve? Matrix inversion has some nasty numerical properties that other methods overcome, and matrix inversion is also slower than most of these methods for some problem. 3. I've run this problem through R so you can see what the inverse of the matrix *should* be. Note that "small" numbers like -4.918979e-17 are zero to within the precision of the operations performed. Also note that I computed the inverse using a process that is more robust than a typical matrix inverse -- it's called the Moore-Penrose Generalized Inverse and is in the MASS library of R as function "ginv". "%*%" is how you do a matrix multiply in R. > m [,1] [,2] [,3] [1,] 0.4177880 5.420376e-01 7.291423e-01 [2,] 0.0000000 6.123032e-17 1.000000e+00 [3,] -0.5420376 4.177880e-01 -2.558129e-17 > ginv(m) [,1] [,2] [,3] [1,] 8.920393e-01 -0.6504235 -1.157331e+00 [2,] 1.157331e+00 -0.8438587 8.920393e-01 [3,] -1.125944e-16 1.0000000 -3.150951e-16 > m%*%ginv(m) [,1] [,2] [,3] [1,] 1.000000e-00 4.857226e-16 -3.852509e-16 [2,] -4.173067e-17 1.000000e+00 -2.604753e-16 [3,] -4.049495e-17 -4.918979e-17 1.000000e-00 ---------------------------------------------------------------------- You can respond by visiting: http://rubyforge.org/tracker/?func=detail&atid=1698&aid=9360&group_id=426