2

I have the following matrix:

[
  [
    16.0,
    23251143833701.0,
    3.3788480598465756e+25,
    4.9101301394821435e+37
  ], 
  [
    23251143833701.0,
    3.3788480598465756e+25,
    4.9101301394821435e+37,
    7.135383882205603e+49
  ],
  [
    3.3788480598465756e+25,
    4.9101301394821435e+37,
    7.135383882205603e+ 49,
    1.0369114809614645e+62
  ], 
  [
    4.9101301394821435e+37,
    7.135383882205603e+49,
    1.0369114809614645e+62,
    1.5068361241656833e+74
  ]
]

The determinant of this matrix is -1.4536774485912138e+135, and if I call regular?, it returns true. However, when I call the method inverse, Ruby raises the following exception:

matrix.rb:1079:in `block in inverse_from': Not Regular Matrix (ExceptionForMatrix::ErrNotRegular)
  from matrix.rb:1069:in `upto'
  from matrix.rb:1069:in `inverse_from'
  from matrix.rb:1061:in `inverse'

What is the problem here?

sawa
  • 165,429
  • 45
  • 277
  • 381
Truncarlos
  • 175
  • 2
  • 7
  • 1
    Seems to be floating point issue. Try to replace your float values with instances of [`BigDecimal`](http://ruby-doc.org/stdlib-2.3.0/libdoc/bigdecimal/rdoc/BigDecimal.html), i.e. use `BigDecimal.new('16.0')` instead of `16.0`. – Stefan Feb 15 '16 at 13:41
  • The error persists. I also believed that it was a floating point issue, but it would be so in case of near-zero values and near-zero determinant, don't you think? – Truncarlos Feb 15 '16 at 14:42

1 Answers1

6

This matrix is numerically inappropriate and cannot be inverted with a general purpose algorithm using floating point numbers. It is singular to the given machine precision:

A =  1e74 * [ 0, 0, 0, 0;
              0, 0, 0, 0;
              0, 0, 0, 0;
              0, 0, 0, 1.5068 ]

Ruby is very optimistic when calculating the determinant. You cannot be certain if -1.45e+135 is a valid result, or a numeric artefact. I tend to trust Matlab and Octave for numerical computation. They both return det(A) = 0.

The determinant of a 4x4 matrix is the sum of products of a permutation of elements, with 4 elements (factors) per product. In the case of this matrix, the smaller products will have exponents around 4*20 = 80, the bigger ones around 4*60 = 240. To calculate sums of the bigger and smaller numbers without loosing too much precision, you would need a floating point number with a mantissa length of at least 240-80 = 160 digits. Single precision float has 24 bits - about 9 or 10 digits, and double precision has 53 bits - about 16 digits. There are probably smarter ways to calculate, but some 160 digits would be needed to calculate with a straight-forward algorithm, e.g. implementing the Leibniz formula.

Even if the inverse matrix was available and precise, it is very questionable if it would be useful in a real-world application.

bogl
  • 1,864
  • 1
  • 15
  • 32
  • 1
    My matrix is an intermediate result of a polynomial regression calculation, using points timestamp-value. The problem was the size of the timestamps (milliseconds), I just converted and normalized them to lower numbers, and the calculation went fine. Thanks! – Truncarlos Feb 17 '16 at 10:11