いけがみです。matrix.rb をつらつらと眺めました。

From: rubikitch <rubikitch / ruby-lang.org>
Subject: [ruby-math:00284] matrix.rb bug?
Date: Fri, 12 May 2000 08:38:44 +0900

> #!/usr/bin/ruby
> load "/usr/local/lib/ruby/1.4/matrix.rb"
> p Matrix[[1,2,3],[-4,-5,-6]].rank
> 
> このスクリプトを実行すると、↓のエラーが出ました。
(エラーの詳細はるびきちさんのメールを参照して下さい)

はい。これは、
rankメソッドで用いる補助行列 a (Array)を transpose した一方で、
列のサイズと行のサイズは転置をとるまえの行列のそれらを利用しているために
起こったバグです。

問題の箇所はここ:
matrix.rb lines 615-
if column_size > row_size
	a = transpose.to_a   # <--- にもかかわらず行と列のサイズはselfのそれ。
else
	a = to_a
end

たとえば、こんなふうにしてみました。
変数を増やした (a_column_size, a_row_size) のはカッコ悪いかもしれませんが...
rankメソッド中の column_size, row_size を配列 a のそれとして、
a_column_size, a_row_size におきかえています。この名前のつけかたは適当です。
# なにせ、今日は1時間しか寝ていないかつ2日ほど水しか口にしていないので ^^;

# 配布するライブラリのメソッドの中にコメントをいれないのは流儀ですか?
# 追えばわかりましたが、非ゼロのピボットが存在しないかもしれないのに
# exists=trueになってる間があるのはちょっと気持ちが悪かったです。
# これが番兵というものでしょうか。

ところで、Matrix.rb の Gauss の消去法は非常に素朴なものですね。
たとえば、浮動小数点の丸め誤差をできるだけ押えるために、
ピボット (rankメソッド中の `q') はできるだけ小さな値にするとか、
列の交換を成分ごとに行なうのではなく、もっと効率の良い方法はないのか、とか。
# たとえば行列が C言語のポインタでされてたら、
# ポインタのつけかえだけで交換が済みますよね。

その意味で、改良の余地は十分にあると思います。

まだ調査をきちんとしてませんが、僕は
・普遍的かつもっとも効率の良い rank を求めるアルゴリズムはなんだろうか
・割算が高価な場合を仮定したときにもっとも効率の良い rank を(以下略)
に興味があります。ガウスの消去法にしても、ピボットの選び方や前処理で
スピードがかなり変わるはずです。
(大きなサイズの行列にならないとその差は有意に出ないでしょうけれど)
それから、割算が高価というのは、
有限体では a/b の計算にいろいろと時間がかかるからです。
体係数有理多項式でも割算が高価なのではないかなと思います。
--
池上 大介
Daisuke IKEGAMI <daisu-ik / is.aist-nara.ac.jp>
奈良先端科学技術大学院大学 情報科学研究科
情報処理学専攻 情報基礎学講座 関研究室

もしかしたら、patchがきちんとあたらないかもしれません。
いつ matrix.rb をいじったか覚えていないのです。
# 原本は保存すべきだったなあ。
やっていることは、
(1) rankメソッドの最初に、 column_size, row_size を
    a_column_size, a_row_size と名前をつけて保存する。
(2) 以下、元の matrix.rb の rankメソッドで、
    column_size, row_size を参照するときは、
	代わりに a_column_size,	a_row_size を参照する。
(3) a = transpose.to_a するときは、 a_column_size, a_row_size も転置する。
の3点のみです。
もし問題があるようでしたら、ruby_srcのtar-ballをとってきて作り直します。

613c613
< 	
---
> 
615c615,618
<     if column_size > row_size
---
> 	a_column_size = column_size
> 	a_row_size = row_size
> 
>     if a_column_size > a_row_size
616a620
> 	  a_column_size, a_row_size = a_row_size, a_column_size
620c624
<     rank = 0
---
> 	rank = 0
624,653c628,657
< 	i = k
< 	exists = true
< 	begin
< 	  if (i += 1) > column_size - 1
< 	    exists = false
< 	    break
< 	  end
< 	end while a[i][k] == 0
< 	if exists
< 	  a[i], a[k] = a[k], a[i]
< 	  akk = a[k][k]
< 	else
< 	  i = k
< 	  exists = true
< 	  begin
< 	    if (i += 1) > row_size - 1
< 	      exists = false
< 	      break
< 	    end
< 	  end while a[k][i] == 0
< 	  if exists
< 	    k.upto(column_size - 1) do
< 	      |j|
< 	      a[j][k], a[j][i] = a[j][i], a[j][k]
< 	    end
< 	    akk = a[k][k]
< 	  else
< 	    next
< 	  end
< 	end
---
> 		i = k
> 		exists = true
> 		begin
> 		  if (i += 1) > a_column_size - 1
> 			exists = false
> 			break
> 		  end
> 		end while a[i][k] == 0
> 		if exists
> 		  a[i], a[k] = a[k], a[i]
> 		  akk = a[k][k]
> 		else
> 		  i = k
> 		  exists = true
> 		  begin
> 			if (i += 1) > a_row_size - 1
> 			  exists = false
> 			  break
> 			end
> 		  end while a[k][i] == 0
> 		  if exists
> 			k.upto(a_column_size - 1) do
> 			  |j|
> 			  a[j][k], a[j][i] = a[j][i], a[j][k]
> 			end
> 			akk = a[k][k]
> 		  else
> 			next
> 		  end
> 		end
655,661c659,665
<       (k + 1).upto(row_size - 1) do
< 	|i|
< 	q = a[i][k] / akk
< 	(k + 1).upto(column_size - 1) do
< 	  |j|
< 	  a[i][j] -= a[k][j] * q
< 	end
---
>       (k + 1).upto(a_row_size - 1) do
> 		|i|
> 		q = a[i][k] / akk
> 		(k + 1).upto(a_column_size - 1) do
> 		  |j|
> 		  a[i][j] -= a[k][j] * q
> 		end
664c668
<     end while (k += 1) <= column_size - 1
---
>     end while (k += 1) <= a_column_size - 1