Matrix bug

An old bug from ruby 1.6 is still unfixed in the ruby 1.8 preview.

irb(main):001:0> require 'matrix’
true
irb(main):002:0> Matrix[[2, 7], [3, 4]].det
-6
irb(main):003:0> require 'mathn’
true
irb(main):004:0> Matrix[[2, 7], [3, 4]].det
-13
irb(main):005:0>

matrix.rb needs mathn to work correctly. The problem is that the algorithm
for the determinant uses Gaussian elimination and assumes ‘/’ has its
mathematical meaning rather than its C-like meaning (the mathematical meaning
is provided by mathn). I see two simple fixes. One is to always require
mathn.rb when matrix.rb is required. The other is to provide an alternative
algorithm for the determinant when a matrix contains Integer entries. The
first is simpler and cleaner, but might break some code.

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

oinkoink+unet@rexx.com (Bret Jolly) wrote in message news:7e7131a1.0304201226.4808d97e@posting.google.com

An old bug from ruby 1.6 is still unfixed in the ruby 1.8 preview.
[…]
This bug affects more than just the determinant (though I think the
underlying problem is the same). Other methods need the mathn library
to work correctly. For example:
irb(main):001:0> require ‘matrix’
true
irb(main):002:0> m = Matrix[[1, 3], [2, 4]]
Matrix[[1, 3], [2, 4]]
irb(main):003:0> n = m.inv
Matrix[[-3, 2], [1, -1]]
irb(main):004:0> mn
Matrix[[0, -1], [-2, 0]]
irb(main):005:0> require ‘mathn’
true
irb(main):006:0> n = m.inv
Matrix[[-2, 3/2], [1, -1/2]]
irb(main):007:0> m
n
Matrix[[1, 0], [0, 1]]

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

irb(main):001:0> require ‘matrix’
=> true
irb(main):002:0> m=Matrix[[1.0,3.0],[2.0,4.0]]
=> Matrix[[1.0, 3.0], [2.0, 4.0]]
irb(main):003:0> n = m.inv
=> Matrix[[-2.0, 1.5], [1.0, -0.5]]
irb(main):004:0> m * n
=> Matrix[[1.0, 0.0], [0.0, 1.0]]

irb(main):001:0> require ‘matrix’
=> true
irb(main):002:0> require ‘rational’
=> true
irb(main):003:0>
irb(main):004:0>
m=Matrix[[Rational.new(1,1),Rational.new(3,1)],[Rational.new(2,1),Rational.new(4,1)]]
=> Matrix[[Rational(1, 1), Rational(3, 1)],
[Rational(2, 1), Rational(4, 1)]]
irb(main):005:0> n = m.inv
=> Matrix[[Rational(-2, 1), Rational(3, 2)],
[Rational(1, 1), Rational(-1, 2)]]
irb(main):006:0> m *n
=> Matrix[[Rational(1, 1), Rational(0, 1)],
[Rational(0, 1), Rational(1, 1)]]

···

On Tuesday, 22 April 2003 at 3:49:23 +0900, Bret Jolly wrote:

oinkoink+unet@rexx.com (Bret Jolly) wrote in message news:7e7131a1.0304201226.4808d97e@posting.google.com

An old bug from ruby 1.6 is still unfixed in the ruby 1.8 preview.
[…]
This bug affects more than just the determinant (though I think the
underlying problem is the same). Other methods need the mathn library
to work correctly. For example:
irb(main):001:0> require ‘matrix’
true
irb(main):002:0> m = Matrix[[1, 3], [2, 4]]
Matrix[[1, 3], [2, 4]]
irb(main):003:0> n = m.inv
Matrix[[-3, 2], [1, -1]]
irb(main):004:0> mn
Matrix[[0, -1], [-2, 0]]
irb(main):005:0> require ‘mathn’
true
irb(main):006:0> n = m.inv
Matrix[[-2, 3/2], [1, -1/2]]
irb(main):007:0> m
n
Matrix[[1, 0], [0, 1]]


Jim Freeze

Too much of a good thing is WONDERFUL.
– Mae West

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> mm.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

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?

[explanation snipped]

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. :slight_smile:

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

My complaint is that I cannot inherit from Matrix
or Vector very easily since they always want to return
themselves and not self.

Should we just write our own and put it in the RAA?
(And add a C version while we’re at it for speed?)

Happy Rubying! Regards, Bret

Yes.

···

On Tuesday, 22 April 2003 at 13:52:07 +0900, Bret Jolly wrote:


Jim Freeze

Heavy, adj.:
Seduced by the chocolate side of the force.

Bret Jolly wrote:

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.)

Ahem, it’s not. The program just “interprets” these numbers in a
different way than you do.
You think of them as real (in the methematical sense), which can be OK
because natural numbers can be trivially incorporated into real numbers.

The progam just uses what you give it: Natural numbers. And consequently
it uses the appropriate operators / (and % probably). Than is, it uses
the (in computer jargon) integer division and modulo operator (often
given as %).

Now speaking in integers 5 / 9 IS 0, while 5 % 9 IS 5.

There’s nothing more to it.

I personally think that you should be as explicit as possible (which is
why I dislike impicit arguments as used in, eg. Perl - note that I’m not
against default values for arguments, though. They are a completely
different thing - and they’re given exlicitly in the definition).
In this case that means that you should probably feed the matrices with
floating point numbers rather than integers.

My (floating point) €0,02

Cheers

Stephan

Now speaking in integers 5 / 9 IS 0, while 5 % 9 IS 5.

There’s nothing more to it.

···

----- Original Message -----
From: “Stephan Kämper” Stephan.Kaemper@Schleswig-Holstein.de


I think it’s clear that the original poster did not need the mysteries of
integer arithmetic explained. :slight_smile:

Now an explanation for you:

Getting back to an earlier post in this thread, we have this:

irb(main):002:0> m = Matrix[[1, 3], [2, 4]]
Matrix[[1, 3], [2, 4]]
irb(main):003:0> n = m.inv
Matrix[[-3, 2], [1, -1]]

Now, however you may feel about integer division, this is just WRONG. That
is not the inverse of that matrix an any sense of the word. That is, in
fact, a totally meaningless answer.

The original poster was not saying that Ruby is buggy. The original poster
was certainly not saying anything about integer arithmetic being buggy.
The original poster was saying that the above is wrong in every situation.
Until you can find any way in which the above is the desired behavior, or is
even meaningful at all, I don’t see why you would feel the need to disagree
with this.


I personally think that you should be as explicit as possible

In this case that means that you should probably feed the matrices with
floating point numbers rather than integers.

?!? Because they are more “explicit”?? You seem to know a lot more about
integer arithmetic than about floating-point arithmetic.

Anyone doing serious floating point computation has to be concerned with
error accumulation, rounding errors, numbers too big or small to fit into a
float… it’s all about minimizing error!

Anyway, here’s another possible solution (instead of just requiring mathn):

It could check when it is loaded if 5/9 is zero or not. If so, it could
generate an error or spit out a warning or something. This would allow one
to use matrix with a different math package, if you wanted. The warning
could even suggest possible solutions (like including mathn instead).

Chris

Jim Freeze jim@freeze.org wrote in message news:20030422055709.A40269@freeze.org

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. :slight_smile:

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/