From: Zinsser <zinsser@...>
Date: 2006-02-02T06:51:53+09:00
Subject: Re: Is this a floating point precision problem?

> Could this be architecture related?  I've tried this both on a pc and 
> a mac and get the same result.  What are you running on?

Hi Todd,

I looked into matrix.rb and found the reason for your problem.
BTW, I also get the wrong result when calling a.inv:

Matrix[[-8.0, 0.0, -0.998844486916645], [0.0, 1.0625, 0.0480592442327619],
[-0.974313737622412, -0.225194894880513, -2.52327742604428e-17]]

It is a floating point precision problem in matrix.rb, or, 
to be more precise, the problem is that the algorithm 
for inverting the matrix is not very sophisticated.


def inverse_from(src)
    size = row_size - 1
    a = src.to_a
    
    for k in 0..size
      if (akk = a[k][k]) == 0
        i = k
        begin
          Matrix.Raise ErrNotRegular if (i += 1) > size
        end while a[i][k] == 0
        a[i], a[k] = a[k], a[i]
        @rows[i], @rows[k] = @rows[k], @rows[i]
        akk = a[k][k]
      end
....

For your matrix, a[1][1] evaluates to -2.33146835171283e-15,
which is close to zero, but does not get caught by 
"if (akk = a[k][k]) == 0". If you replace this version with 

def inverse_from(src)
    size = row_size - 1
    a = src.to_a

    puts a
    
    for k in 0..size
      if (akk = a[k][k]) > -1.0e-12 && akk < 1.0e-12
        i = k
        begin
          Matrix.Raise ErrNotRegular if (i += 1) > size
        end while a[i][k] == 0
        a[i], a[k] = a[k], a[i]
        @rows[i], @rows[k] = @rows[k], @rows[i]
        akk = a[k][k]
      end
....

you get the correct result. Perhaps someone could implement 
this algorithm with a better pivoting strategy?

Your matrix seems to be orthonormal, aka a 3-D rotation matrix.
For this kind of matrix, the inverse is identical to the transpose,
which is also much faster to compute.

Timo Zinsser