[SciPy-user] determinants
David Goldsmith
David.L.Goldsmith@noaa....
Thu Sep 13 16:55:38 CDT 2007
Thanks!
DG
Robert Kern wrote:
> 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.
>
>
More information about the SciPy-user
mailing list