Bugs item #9360, was opened at 2007-03-18 14:41
You can respond by visiting: 
http://rubyforge.org/tracker/?func=detail&atid=1698&aid=9360&group_id=426

Category: Standard Library
Group: 1.8.6
Status: Open
Resolution: None
Priority: 3
Submitted By: Stefan Mühlebach (muehlebach)
Assigned to: Nobody (None)
Summary: Matrix inverse problem

Initial Comment:
There is a problem with the Matrix.inverse method.
Given the following matrix

    m = Matrix[[0.417787968829298, 0.542037605627772, 0.729142268138943],
    [0.0, 6.12303176911189e-17, 1.0],
    [-0.542037605627772, 0.417787968829298, -2.55812900589452e-17]]

the result of

    m.inverse * m

is

    [[1.0, 0.0, 0.650423523448542],
    [5.55111512312578e-17, 1.0, 0.84385869292008],
    [0.0, 9.24446373305873e-33, 1.0]]

but it should be

    [[1.0, 0.0, 0.0],
    [0.0, 1.0, 0.0],
    [0.0, 0.0, 1.0]]

or close to this!


----------------------------------------------------------------------

Comment By: Ed Borasky (znmeb)
Date: 2007-03-18 15:27

Message:
Questions:

1. Have you loaded "mathn"? If not, load it and try again. The "Matrix" routines are more accurate with "mathn" loaded.

2. What is the problem you're actually trying to solve? Matrix inversion has some nasty numerical properties that other methods overcome, and matrix inversion is also slower than most of these methods for some problem.

3. I've run this problem through R so you can see what the inverse of the matrix *should* be. Note that "small" numbers like -4.918979e-17 are zero to within the precision of the operations performed. Also note that I computed the inverse using a process that is more robust than a typical matrix inverse -- it's called the Moore-Penrose Generalized Inverse and is in the MASS library of R as function "ginv". "%*%" is how you do a matrix multiply in R.

> m
           [,1]         [,2]          [,3]
[1,]  0.4177880 5.420376e-01  7.291423e-01
[2,]  0.0000000 6.123032e-17  1.000000e+00
[3,] -0.5420376 4.177880e-01 -2.558129e-17
> ginv(m)
              [,1]       [,2]          [,3]
[1,]  8.920393e-01 -0.6504235 -1.157331e+00
[2,]  1.157331e+00 -0.8438587  8.920393e-01
[3,] -1.125944e-16  1.0000000 -3.150951e-16
> m%*%ginv(m)
              [,1]          [,2]          [,3]
[1,]  1.000000e-00  4.857226e-16 -3.852509e-16
[2,] -4.173067e-17  1.000000e+00 -2.604753e-16
[3,] -4.049495e-17 -4.918979e-17  1.000000e-00



----------------------------------------------------------------------

You can respond by visiting: 
http://rubyforge.org/tracker/?func=detail&atid=1698&aid=9360&group_id=426