[ruby-dev:27730] Re: Matrix class is broken without mathn

原です。

ちょっと間が空いてしまいました。すいません。Integer成分の
determinantの話の続きです。

In [ruby-dev:27700]
>けいじゅ＠いしつかです.

>>でも、quoを使った方が現状よりは（クレームの頻度の点で）かなりまし
>>じゃないですかね。実際にdetを計算したいのはFloatの行列が多いらしい
>>し。
>
>ましなのは理解しますが, Rationalをrequireするよりましな理由がわかりま
>せんが?

私が誤解しているのか、どうも石塚さんの話がぴんとこないのです
が、、、私は、現行ではRationalをrequireするだけでなく、更に
各成分をRational化してdetに放り込まなくてはならないのが問題
だと思うのです。

>det計算中Rationalをrequireする場合,
>すべてがIntegerの場合, 結果は必ずInitegerになり,
>行列中1要素でもFloatが入っている場合, 結果は必ずFloatになり, 
>Rationalになることはありません.

今回問題になっているのは、整数行列の行列式はどうあつかうべき
か、ですよね？

>ですので, (内部の計算がどうなっているかは別として), それぞれの
>クラスに応じた自然な結果が出ることが必ず保証できます.

細かい話ですが、整数行列の成分をRational(x, 1)で全て有理数
化してdetを計算させると、結果はRationalですよね。mathn.rbを
導入するとIntegerですけど。

>あと, パフォーマンスのことを気にするなら, 最初っから要素は
>Floatにすべきだと思います.

私がquoを使った方がいいと思うのは、整数行列のdetを計算したと
きとんでもない値が出るのを防げるからという理由で、ただその一
点に尽きます。

それに対して石塚さんはrational.rbをrequireして各成分、
Rational()でラップしてくれ、と主張しているわけですよね。
私はそれは面倒すぎるのではないかと、そして今後も忘れて
意図しない結果に悩む人が大勢出るのではないかと思います。
それならちょっとかっこ悪くても「quoでFloat」でいいのでは
ないでしょうか。

現行をquo式にしたメリットデメリットまとめると、

--------------------------------------------------------------
成分               | 現行               | quo式               
-------------------+--------------------+---------------------
全てInteger        | NG                 | 戻値:Float          
                   |                    | 誤差が気になることも
-------------------+--------------------+---------------------
全てInteger        | NG                 | 戻値:Rational       
rationa.rb導入のみ |                    | 遅い                
-------------------+--------------------+---------------------
全てRational       | 戻値:Rational      | 戻値:Rational       
                   | 遅さは気にならない | 遅さは気にならない  
-------------------+--------------------+---------------------
全てFloat          | 戻値:Float         | 戻値:Float          
                   | 誤差は気にならない | 誤差は気にならない  
-------------------+--------------------+---------------------

となります。「遅さは気にならない」とは、Rational計算である
のが明白なので遅いことは覚悟が出来ている、、という意味です。
:-) 
quo式は、rational.rbをrequireしてるか否かで結果が違うのが確
かに気持ちが悪いのですが、それ以上にとにかく現行の"NG"は良く
ないのではないか。

quo式で困るとしたら、rational.rbをrequireしている状態で、整
数行列のdetをFloatで高速に欲しい人（いるのか？）ですが、その
人には各成分Float()してからquo式detを計算してもらうことにし
ます。

どうでしょう。


更にもう一つ提案なんですが、整数行列のdetの計算には、互除法
を使ったそこそこ速いのがあります：

  def determinant_e
    size = row_size - 1
    a = to_a

    det = 1
    0.upto(size) do |k|
      (k + 1).upto(size) do |i|
        until a[i][k].zero?
          q = a[k][k] / a[i][k]
          k.upto(size) do |j|
            a[k][j] -= a[i][j] * q
          end
          a[k], a[i] = a[i], a[k]
          det *= -1
        end
      end
      if a[k][k].zero?
        return 0
      else
        det *= a[k][k]
      end
    end
    det
  end

これは、整数行列をRatinal()してから、現行detを計算するより
かなり速いです。次は15次行列の整数成分の行列式を5回計算させ
たのにかかった時間です。

determinant         : 4.278sec.
determinant_e       : 0.329sec.

上段が各成分をRational()して現行のdetによる計算です。
拡張ライブラリのrational.soを使っても1.4秒でした。

このdeterminant_eの面白いのは、IntegerとRationalの混合した
成分を持つ行列でも正しく動作する所です。まとめると、

------------------------------------------------------------
成分               | quo式                | determinant_e
-------------------+----------------------+-----------------
全てInteger        | 戻値:Float           | 戻値:Integer    
                   | 誤差が気になることも | やや遅い        
-------------------+----------------------+-----------------
全てInteger        | 戻値:Rational        | 戻値:Integer    
rationa.rb導入のみ | 遅い                 | やや遅い        
-------------------+----------------------+-----------------
全てRational       | 戻値:Rational        | 戻値:Rational   
                   | 気にならない遅さ     | 気にならない遅さ
-------------------+----------------------+-----------------
全てFloat          | 戻値:Float           | NG              
                   | 気にならない誤差     |                 
-------------------+----------------------+-----------------

そこで提案なんですが、Matrix のdetはquo式にし、と同時に
determinant_e（det_iという名前もいいですが）も定義しておいて、

  【整数行列を扱う場合の注意】
    整数行列の行列式を計算させるときは、Matrix#det_iを使いましょ
    う。Matrix#detは、rational.rbをrequireしていない状態でFloat、
    requireしている状態でRationalを返します。

と、マニュアルに記述するのです。

どうでしょう？

ちなみに、determinant_eは[ruby-math:625]のコードをシンプルに
したものです。他にも幾つかバリエーションを試してみたのですが、
結局[ruby-math:625]のコードが速いのでこちらの方がいいかもし
れません。

テストスクリプトを載せます。

#!/usr/local/bin/ruby
require "matrix"
require "rational"

class Matrix
  def determinant_e
    size = row_size - 1
    a = to_a

    det = 1
    0.upto(size) do |k|
      (k + 1).upto(size) do |i|
        until a[i][k].zero?
          q = a[k][k] / a[i][k]
          k.upto(size) do |j|
            a[k][j] -= a[i][j] * q
          end
          a[k], a[i] = a[i], a[k]
          det *= -1
        end
      end
      if a[k][k].zero?
        return 0
      else
        det *= a[k][k]
      end
    end
    det
  end

  alias det_e determinant_e

  #[ruby-math:625]
  def determinant_625
    return 0 unless square?
    
    size = row_size - 1
    a = to_a
    
    det = 1
    k = 0
    begin 
      if a[k][k].zero?
        i = k
        begin
          return 0 if (i += 1) > size
        end while a[i][k].zero?
        a[i], a[k] = a[k], a[i]
        det *= -1
      end
      (k + 1).upto(size) do |i|
        q = a[i][k] / a[k][k]
        k.upto(size) do |j|
          a[i][j] -= a[k][j] * q
        end
        unless a[i][k].zero?
          a[i], a[k] = a[k], a[i]
          det *= -1
          redo
        end
      end
      det *= a[k][k]
    end while (k += 1) <= size
    det
  end
end

class Test
  attr_reader :value
  def initialize(name)
    @name = name
    @value = nil
    @time = 0.0
  end

  def stime(v = nil)
    t = Time.now
    @value = yield
    @time += Time.now - t
    if v && v != @value
      puts "Test Failed: #{@value} != #{v}"
      exit
    end
  end

  def to_s
    sprintf("%-30s: %3.3fsec.", @name, @time)
  end
end

def test(dim, num, range = 100)
  tr = Test.new("determinant(rational)")
  te = Test.new("determinant_e(int)")
  t1 = Test.new("determinant_e(int & rational)")
  tm = Test.new("determinant_625(int)")
  tf = Test.new("determinant(float)")

  num.times do
    a = (0...dim).map{(0...dim).map{rand(2*range+1)-range}}
    m = Matrix[*a]
    mr = Matrix[*a.map{|b| b.map{|c| Rational(c, 1)}}]
    m1 = Matrix[*a.map{|b| b.map{|c|
          rand(2).zero? ? Rational(c, 1) : c}}]
    mf = Matrix[*a.map{|b| b.map{|c| Float(c)}}]

        
    tr.stime{Matrix[*mr].det}
    te.stime(tr.value){Matrix[*m].determinant_e}
    t1.stime(tr.value){Matrix[*m1].determinant_e}
    tm.stime(tr.value){Matrix[*m].determinant_625}
    tf.stime(){Matrix[*mf].determinant}
  end

  puts "--- #{dim}-dim, #{num}times ---"
  puts tr, te, t1, tm, tf
  puts
end

test(5, 500)
test(10, 20)
test(15, 5)
test(20, 2)


以下、テストの結果：

--- 5-dim, 500times ---
determinant(rational)         : 4.503sec.
determinant(rational)         : 0.647sec. (with rational.so)
determinant_e(int)            : 1.158sec.
determinant_e(int & rational) : 7.313sec.
determinant_625(int)          : 0.703sec.
determinant(float)            : 0.193sec.

--- 10-dim, 20times ---
determinant(rational)         : 3.058sec.
determinant(rational)         : 0.786sec. (with rational.so)
determinant_e(int)            : 0.421sec.
determinant_e(int & rational) : 4.223sec.
determinant_625(int)          : 0.251sec.
determinant(float)            : 0.077sec.

--- 15-dim, 5times ---
determinant(rational)         : 4.278sec.
determinant(rational)         : 1.390sec. (with rational.so)
determinant_e(int)            : 0.329sec.
determinant_e(int & rational) : 5.347sec.
determinant_625(int)          : 0.238sec.
determinant(float)            : 0.027sec.

--- 20-dim, 2times ---
determinant(rational)         : 5.384sec.
determinant(rational)         : 2.208sec. (with rational.so)
determinant_e(int)            : 0.884sec.
determinant_e(int & rational) : 6.511sec.
determinant_625(int)          : 0.766sec.
determinant(float)            : 0.024sec.