[ ruby-Bugs-9360 ] Matrix inverse problem

Bugs item #9360, was opened at 2007-03-18 23: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: Peter Vanbroekhoven (calamitas)
Date: 2007-03-19 03: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-19 03: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-19 02: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-19 01: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-19 01: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-19 01: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-19 00: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-19 00: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