From: oinkoink+unet@... (Bret Jolly)
Date: 2003-04-22T13:52:07+09:00
Subject: Re: Matrix bug

Jim Freeze <jim@freeze.org> wrote in message news:<20030421154222.A38737@freeze.org>...

> Is the problem that you are asking for 5/9 not to be zero?

   The problem is that, when mathn is not required, matrix.rb
allows 5/9 to be zero and this is just wrong. (Yes, I know C;
it's my job, but this is still mathematically wrong.) It leads 
to the WRONG results, as I've shown, for determinants and 
inverses and undoubtably for nearly everything else.  This is
a big bug. It can be fixed by forcing mathn to be required
when matrix is required (perhaps by incorporating matrix in 
mathn.rb -- currently mathn.rb requires matrix).

   The determinant could be fixed by changing the algorithm,
but the inverse is not so easy.  The inverse of an integer
matrix can have non-integer entries.  You could replace the
integer entries by floats, but this would be an abomination.
Exact operations on exact data should lead to exact results.
Theory yields exact matrices that are non-singular but very
ill-conditioned. (Perhaps the best known is the Hilbert matrix,
which arises naturally in fitting least-square polynomials to
data.  The Hilbert matrix has rational entries, but its
inverse has integer entries.) Floating point arithmetic
on such beasts will yield wildly wrong results, no matter how
careful you are with your numerics.

   If you are curious about these matrices, here is some code.
Notice that I require mathn!

# hilbert.rb
# Hilbert matrix of order n
require 'mathn'

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

irb(main):001:0> require 'hilbert'
true
irb(main):002:0> m = hilbert(5)
Matrix[[1, 1/2, 1/3, 1/4, 1/5], [1/2, 1/3, 1/4, 1/5, 1/6],
 [1/3, 1/4, 1/5, 1/6, 1/7], [1/4, 1/5, 1/6, 1/7, 1/8],
 [1/5, 1/6, 1/7, 1/8, 1/9]]
irb(main):003:0> m.det
1/266716800000  # damn near zero => matrix is nearly singular
irb(main):004:0> m.inv
Matrix[[25, -300, 1050, -1400, 630], [-300, 4800, -18900, 26880,
-12600], [1050, -18900, 79380, -117600, 56700], [-1400, 26880,
-117600, 179200, -88200], [630, -12600, 56700, -88200, 44100]]
# NB: Integer entries!

  Currently, matrix.rb does not work properly with integer 
matrices unless mathn is included. Notice that you have to 
require mathn; it is not enough just to require rational.rb:
irb(main):001:0> require 'matrix'
true
irb(main):002:0> require 'rational'
true
irb(main):003:0> m = Matrix[[1, 2], [3, 4]]
Matrix[[1, 2], [3, 4]]
irb(main):004:0> m*m.inv
Matrix[[0, -1], [-2, -1]]
irb(main):005:0> require 'mathn'
true
irb(main):006:0> m*m.inv
Matrix[[1, 0], [0, 1]]

Happy Rubying! Regards, Bret