[SciPy-user] determinants

Robert Kern robert.kern@gmail....
Thu Sep 13 16:45:35 CDT 2007


David Goldsmith wrote:
> Anne Archibald wrote:
>> On 13/09/2007, David Goldsmith <David.L.Goldsmith@noaa.gov> wrote:
>>   
>>> Anne Archibald wrote:
>>>     
>>>> If you have problems with determinants becoming excessively large, you
>>>> may be able to circumvent them by computing the log of the
>>>> determinant. The easiest way to do this is to use LU decomposition:
>>>>
>>>> P,L,U = scipy.linalg.lu(M)
>>>> d = sum(log(abs(diag(L))))
>>>>
>>>> Of course you lose track of the sign doing this (P may be either an
>>>> even or odd permutation, though det should be reliable and efficient
>>>> on it).
>>>>       
>>> Of course, one can keep track of the sign by cumprod(sgn(diag(L))),
>>> yes?  (Sorry, I don't know the numpy functions for these off hand, but I
>>> assume they exist, yes?)
>>>     
>> In my first draft I suggested this (prod(sign(diag(L)))), but
>> unfortunately P may contribute a factor of -1, so you have to extract
>> its determinant as well (unless there's some more clever way to get
>> the sign of a permutation matrix? some sum involving the positions of
>> the ones modulo 2 ought to do it, but it's been a while since I did
>> this kind of combinatorics).
>>   
> Right, (clearly) same here. ;-)
> 
> However, this suggests something that maybe should be implemented inside 
> det's "black box"?  (Obviously, if det is being used inside a formula, 
> the function can't simply return log(det) w/out some manner of user 
> notification.  Perhaps a custom exception and/or a second function, e.g. 
> "logdet", to which the user is referred if abs(det) returns inf?  Just a 
> suggestion.)

I would just implement logdet() separately and not try to switch between the two
automatically. It should be fairly straightforward. You can look in
scipy/linalg/src/det.f and scipy/linalg/basic.py for the relevant pieces of code.

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


More information about the SciPy-user mailing list