From: noreply@... Date: 2007-03-20T12:11:41+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: Closed 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-19 19:11 Message: Sure, I'll come up with some. I can generate/find some pretty nasty matrices that will mess up almost any inverse calculation and still technically have full rank. They should work if you define them as rationals rather than doubles, though, and bring in the library using "require mathn". By the way, I've proposed adding some more algorithms (LU, QR and singular value decompositions, eigensolutions, etc.) as a Summer of Code project. ---------------------------------------------------------------------- Comment By: Akinori MUSHA (knu) Date: 2007-03-19 18:50 Message: Committed to trunk and ruby_1_8, thanks! Now, does anybody care to write test cases? ---------------------------------------------------------------------- Comment By: Peter Vanbroekhoven (calamitas) Date: 2007-03-18 18:33 Message: Sorry for replying to myself *again*. Next time I should first test a patch before submitting. Fixed version of second patch below. (Note the added akk = a[k][k].) Peter --- lib/matrix.rb.old 2007-03-19 01:46:44.000000000 +0100 +++ lib/matrix.rb 2007-03-19 03:26:38.000000000 +0100 @@ -590,15 +590,21 @@ 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 + i = k + akk = a[k][k].abs + for j in (k+1)..size + v = a[j][k].abs + if v > akk + i = j + akk = v + end + end + Matrix.Raise ErrNotRegular if akk == 0 + if i != k a[i], a[k] = a[k], a[i] @rows[i], @rows[k] = @rows[k], @rows[i] - akk = a[k][k] end + akk = a[k][k] for i in 0 .. size next if i == k ---------------------------------------------------------------------- Comment By: Peter Vanbroekhoven (calamitas) Date: 2007-03-18 18:23 Message: This version seems to be slightly faster on the OP's example (9% on my machine). In my opinion the version with the sort_by is more elegant though. If only 1.8 had 1.9's max_by... Peter --- lib/matrix.rb.old 2007-03-19 01:46:44.000000000 +0100 +++ lib/matrix.rb 2007-03-19 02:56:38.000000000 +0100 @@ -590,14 +590,19 @@ 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 + i = k + akk = a[k][k].abs + for j in (k+1)..size + v = a[j][k].abs + if v > akk + i = j + akk = v + end + end + Matrix.Raise ErrNotRegular if akk == 0 + if i != k a[i], a[k] = a[k], a[i] @rows[i], @rows[k] = @rows[k], @rows[i] - akk = a[k][k] end for i in 0 .. size ---------------------------------------------------------------------- Comment By: Peter Vanbroekhoven (calamitas) Date: 2007-03-18 17:21 Message: The patch below should "fix" the inverse method. It adds partial pivoting to the Gauss-Jordan algorithm, making it stable. (I used Ruby 1.8.5 (2006-08-22).) I don't think there is an error in the algorithm, and the matrix is not ill-conditioned (a condition number of about 2.5 isn't so bad). Looks to me like BigDecimal uses Gaussian elimination with partial pivoting too, or at least that's what the comments say. It includes some additional scaling, but this may be specific to BigDecimal (precision works differently for BigDecimals and floats). Peter --- lib/matrix.rb.old 2007-03-19 01:46:44.000000000 +0100 +++ lib/matrix.rb 2007-03-19 01:46:06.000000000 +0100 @@ -590,15 +590,13 @@ 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 + i = (k..size).sort_by { |j| a[j][k].abs }[-1] + Matrix.Raise ErrNotRegular if a[i][k] == 0 + if i != k a[i], a[k] = a[k], a[i] @rows[i], @rows[k] = @rows[k], @rows[i] - akk = a[k][k] end + akk = a[k][k] for i in 0 .. size next if i == k ---------------------------------------------------------------------- Comment By: Stefan M��hlebach (muehlebach) Date: 2007-03-18 16:56 Message: What the hell is 'R'? The equation system I'm trying to solve looks like this: 1) n * ray = cos(phi) 2) ray * ray' = cos(2*phi) 3) ray' * (n x ray) = 0 Where phi is a given scalar, n and ray' are given vectors with three elements, '*' is the scalar product and 'x' is the vector product. The elements of ray are to be computed. The equations are part of a nice Java-Applet (http://muehlebach.dyndns.org/~stefan/Egg/EggApplet.html - You need the J3D-Library to run the applet! Use the mouse and the cursor keys.) I just installed GSL and try to get into it... ---------------------------------------------------------------------- Comment By: Ed Borasky (znmeb) Date: 2007-03-18 16:42 Message: OK ... I found the FORTRAN code, which is as difficult to read as what's in "Matrix". :) The algorithms look close; the FORTRAN version claims it is using a Gauss-Jordan algorithm. Bear in mind that the FORTRAN in question is from the mid 1960s, and that the Gauss-Jordan algorithm has not been used for matrix inversion for decades. So there are two possibilities that I can see: your matrix is so poorly conditioned that Gauss-Jordan is failing to invert it, or the programmer of Matrix made a mistake in coding the Gauss-Jordan. I suspect the latter, since some of the rows and columns look right and some of them don't. In either case, my recommendation would be to replace the Gauss-Jordan with something "modern", like an LU decomposition with partial pivoting and multiple right-hand sides. Because of Ruby's "duck typing", you can probably lift the LU decomposition code right out of the solver that's in BigDecimal. :) As soon as someone picks this bug up, I'll be happy to provide references on computational linear algebra. ---------------------------------------------------------------------- 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