Error in linalg.inv ??

[Robert Kern]

On 2009-06-06 17:09, John Machin wrote:
> Robert Kern<robert.kern<at>  gmail.com>  writes:
>
>> On 2009-06-06 00:34, Ajith Kumar wrote:
>
>>> from numpy import *
>> Please ask numpy questions on the numpy mailing list, not here:
>>
>>     http://www.scipy.org/Mailing_Lists
>>
>>> a = arange(1.,10.)
>>> b = reshape(a, [3,3])
>>> c = linalg.inv(b)
>>>
>>> And the result is
>>>
>>> [[ 3.15221191e+15 -6.30442381e+15 3.15221191e+15]
>>> [ -6.30442381e+15 1.26088476e+16 -6.30442381e+15]
>>> [ 3.15221191e+15 -6.30442381e+15 3.15221191e+15]]
>
> (gibberish)
>
>> You have a very singular matrix (2*a[1] - a[0] == a[2]). You cannot invert it
>> numerically and expect sensible results.
>
> Is raising an exception (as documented
> (http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.inv.html)) not
> a "sensible result"?

That's not the way I was using the phrase, but you have a point.

We depend on the underlying LAPACK routine to tell us if the array is singular 
or not. It will inform us if it stopped prematurely because of singularity. In 
this case, I think it happens to complete simply because it is so small and the 
subroutine did not detect the singularity. A (5,5) matrix constructed the same 
way will raise the exception, for example.

However, the condition number of a matrix is not a binary property. A matrix can 
be poorly conditioned but not precisely singular, and the inversion routine will 
give you a correct result (within the limits of floating point arithmetic) 
rather than an exception, but you will still get mostly nonsense out if you try 
to use that result. You still need to apply care in interpreting the result and 
not rely on our throwing an exception even if we could catch all singular inputs.

-- 
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless enigma
  that is made terrible by our own mad attempt to interpret it as though it had
  an underlying truth."
   -- Umberto Eco