Re: Matrix bug

Jim Freeze <[e-mail]> wrote in message news:<[e-mail]>...

> Thanks for the fine explanation. I have been able to
> work around such a 'problem' by forcing everything to 
> a float and have not considered it to be too much of
> an abomination. :)

> How would you solve this? Would you just want Matrix 
> to convert to float internally for you before an inverse?

No! No! No!   Sorry for shouting, but I'm an excitable
piggy.  A matrix times its inverse is the identity matrix:
ones on the main diagonal and zeros elsewhere.  So if
we multiply a matrix times by its inverse and get large
off-diagonal entries, we know that something is severely
wrong.  Consider a floating-point variant of the code
I gave yesterday, with a test for large off-diagonal
entries.

# hilbert_float.rb
# Hilbert matrix of order n, converted to floating point
require 'mathn'

def hilbert_float(n)
  row_array = []
  (1..n).each {|r|
     row = []
     (r..(n+r-1)).each{|s| row.push((1/s).to_f)}
     row_array.push(row)
   }
  Matrix[*row_array]
end

class Matrix
  # find off-diagonal element of maximum magnitude
  def biggest_offdiag
    biggest = 0
    (0...row_size).each {|i|
      (0...column_size).each {|j|
        if i != j and biggest.abs < self[i, j].abs
          biggest = self[i, j]
        end
      }
    }
    biggest
  end
end
# end of file

irb(main):001:0> require 'hilbert_float'
true
irb(main):002:0> ((m = hilbert_float(15))*m.inv).biggest_offdiag
-569.0
irb(main):003:0> ((m = hilbert_float(18))*m.inv).biggest_offdiag
-4056.125
irb(main):004:0> ((m = hilbert_float(20))*m.inv).biggest_offdiag
-14119.25
irb(main):005:0> ((m = hilbert_float(25))*m.inv).biggest_offdiag
-180865.0
irb(main):006:0> ((m = hilbert_float(35))*m.inv).biggest_offdiag
-6862824.5

   Cowabunga!  But the hilbert.rb I wrote with rational entries
(and requiring mathn) gave the exact answer.  Thus converting
to Float is not a solution.

    I'd better explain what mathn does, since English Ruby
documentation pretty much ignores it.  First mathn requires
rational.rb, complex.rb, and matrix.rb.  It then performs some
magic to make these and the built-in Numeric classes play nice
with one another.  Effectively, one can now regard Rational as
a subfield of Complex, and Integer as a subring of Rational,
just as any mathematician would want.  The expression 5/9
now evaluates to Rational(5, 9) rather than 0.   You can use
the methods .div or .divmod for whole-number division.  The
methods in matrix.rb require this behavior in order to work
properly.  C-like integer division with '/' gives the wrong
answers.

Regards, Bret
Bret Jolly (Jou Lii Bair) 
http://www.rexx.com/~oinkoink/