Michael Neumann said:
Ah thanks.... would be great to see this in Ruby too (without breaking
existing code).
It would *not* be great to see this in Ruby: the quotient of two integers
should be a rational, not a float. The behavior given by the mathn library
is the correct one. If the quotient of two rationals is a float, you get
severe errors using the matrix standard library with ill-conditioned matrices.
The Hilbert matrix, which arises in the the theory of least-square fits, is
a good example.
First recall that you *must* require mathn if you use the matrix library.
Here's what happens if you forget:
irb(main):001:0> require 'matrix'
true
irb(main):002:0> m = Matrix[[2, 3], [3, 5]]
Matrix[[2, 3], [3, 5]]
irb(main):003:0> m.det
4 # WRONG
irb(main):004:0> m.inv
Matrix[[0, 0], [-1, 0]] # WRONG
irb(main):005:0> m*m.inv
Matrix[[-3, 0], [-5, 0]] # WRONG: m*m.inv must be the identity matrix
irb(main):006:0> require 'mathn'
true
irb(main):007:0> m.det
1 # correct
irb(main):008:0> m.inv
Matrix[[5, -3], [-3, 2]] # correct
irb(main):009:0> m*m.inv
Matrix[[1, 0], [0, 1]] # correct
Now let's look at what happens if you use float results for integer division.
I define two functions: hilbert and hilbert_float, which create a Hilbert
matrix and its floating-point version, respectively.
# hilbert.rb
# Hilbert matrix of order n
require 'mathn' # ESSENTIAL!
# exact (rational entries)
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
# floating-point entries
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
A matrix times its inverse is the identity matrix -- 1s on the main
diagonal and 0s elsewhere. So for any invertible matrix m,
(m*m.inv).biggest_offdiag should be 0. Let's try it.
irb(main):001:0> require 'mathn'
true
irb(main):002:0> require 'hilbert'
true
irb(main):003:0> ((m=hilbert(15))*m.inv).biggest_offdiag
0
irb(main):004:0> ((m=hilbert_float(15))*m.inv).biggest_offdiag
-569.0 # INCREDIBLY WRONG!
The quotient of two integers is an exact rational. To make the
quotient default to (inexact) floats can lead to appalling numerical
errors as I have just shown.
Regards, Bret
