logdet of symmetric positive definite matrix

[argriffing]

I'm trying to use algopy to compute this, but it only has det but not logdet built in. I see that algopy det is computed by taking a square of product of diagonals of a cholesky decomposition of a positive definite matrix. Maybe the logdet could be analogously computed by taking twice the sum of logs of the diagonal of this decomposition?

I'm not sure how the UTPM math works in this case. By the way, I cannot just use logdet := log(det) because it is less stable because the intermediate determinant might be huge even if logdet is reasonably sized.

Edit: I'm interested in this function because it is used in log likelihoods and Kullback-Leibler divergences for multivariate normal distributions.

[b45ch1]

I'll have to fix the implementation of det: It should be computed by a QR decomposition or a singular value decomposition, not by the Cholesky factorization.

about logdet: If you use the QR decomposition, i.e.

A = QR

then the determinant is

det(A) = prod(diag(R))

and hence

logdet(A) = log(prod(diag(R))) = sum(log(diag(R)))

This method should be "stable" in the sense that no huge intermediates show up. Should we add this function to algopy?

If we found an analytical recurrence for the Taylor propagation of the determinant, i.e.

c_0 + c_1 t + c_2 t^2 + ... = det(A_0 + A_1 t + A_2 t^2 + ....)

we could also use that. But I can only find first-order results in the literature.

EDIT:

it should be det(A) = det(Q) * prod(diag(R)) and logdet(A) = logdet(QS) + logdet(SR), where S = sign(diag(R))

akin to http://www.mathworks.com/matlabcentral/fileexchange/22026-safe-computation-of-logarithm-determinat-of-large-matrix/content/logdet.m

[argriffing]

On one hand, logdet is an extremely important function in statistics because it is used for multivariate normal log likelihood. So for statistics it is almost more important than det. But on the other hand, it is messy for matrices that are not positive definite. Maybe algopy could follow cvxpy which has just added a det_rootn function that assumes a positive definite input matrix, and then logdet is implemented as n * det_rootn(A). If there were not a plain det function and if det_rootn and logdet had big warnings in the docs about the input assumptions, then this might be less confusing.

Maybe det could be implemented with QR allowing unrestricted input, and det_rootn and logdet could be implemented with Cholesky decomposition and assume positive definite input matrices?

[b45ch1]

I missed an important point: I was, wrongly, assuming that 1==det(Q), but this is not guaranteed by the numpy/lapack implementation, only dot(Q.T,Q) = Id and thus det(Q) is potentially -1! So the trick det(A) = det(Q)*det(R) = det(R) does not work.

If we called the underlying Lapack routine directly we could compute det(Q) in an inexpensive manner, see http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=2&t=1741 . Getting it directly from Q seems to be difficult.

In numpy, the determinant is computed via an LU decomposition. I'd suggest to add this factorization to AlgoPy and use it to compute det and logdet.

[b45ch1]

The QR factorization in scipy offers the possibility to compute the QR decomposition in 'raw' mode:

Definition: scipy.linalg.qr(a, overwrite_a=False, lwork=None, mode='full', pivoting=False)
Docstring:

Compute QR decomposition of a matrix.

Calculate the decomposition ``A = Q R`` where Q is unitary/orthogonal
and R upper triangular.

Parameters
----------
a : array, shape (M, N)
    Matrix to be decomposed
overwrite_a : bool, optional
    Whether data in a is overwritten (may improve performance)
lwork : int, optional
    Work array size, lwork >= a.shape[1]. If None or -1, an optimal size
    is computed.
mode : {'full', 'r', 'economic', 'raw'}
    Determines what information is to be returned: either both Q and R
    ('full', default), only R ('r') or both Q and R but computed in
    economy-size ('economic', see Notes). The final option 'raw'
    (added in Scipy 0.11) makes the function return two matrixes
    (Q, TAU) in the internal format used by LAPACK.

I can't find a function in scipy that converts (Q,TAU) to Q and R, though.

[b45ch1]

The code around https://github.com/scipy/scipy/blob/master/scipy/linalg/decomp_qr.py#L163 seems to be doing what we need.

[argriffing]

Thanks for adding these improvements to algopy! Personally I only use logdet and only for positive definite matrices. But I understand wanting to have a simple interface that doesn't require the user to say anything about the definiteness of their matrices, and to have det work correctly if you are keeping it. I don't know much about det calculation for matrices that are not positive definite.

[b45ch1]

I have added algopy.det, algopy.logdet, algopy.lu, algopy.lu2, algopy.prod and a few supplementary functions to algopy. algopy.det and algopy.logdet both use algopy.lu2 internally.

There are corresponding unit tests, but I did not have the time for thorough testing. Could you confirm that their implementations work correctly?

[b45ch1]

Closing due to inactivity.