まつもと ゆきひろです

rubyforgeバグレポートに以下のようなレポートが来てますが、私
には判断できません。どうでしょう。石塚さん?

# コメントは編集してあります。パッチは(3つある)最終版のもの
# です。

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

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