topic Re: EM-algorithm in SAS/IML Software and Matrix Computations

EM-algorithm

renc — Wed, 24 Oct 2012 13:58:15 GMT

I did an EM-algorithm for latent variable model with missing biomarkers and coviates in SAS IML. There is an overflow error in the DET of a 159X159 matrix. I have to use Log-likelihood function(LF) to control the iteration times and it includes the determinant.

What are the possible approaches to solve this problem?

Thanks for any help in advance!

Re: EM-algorithm

Rick_SAS — Wed, 24 Oct 2012 14:18:24 GMT

The determinant is a polynomial in the entries, and is therefore known to be unstable for large matrices.

Are you using the determinant to find out if a matrix is singular? If so, there are better (more stable) ways to determine singularity.

Or does the determinant appear in the LL formula and you need to compute it explicitly? If so, compute the log-determinant (which should be stable) and exponentiate, assuming the log-det is less than CONSTANT('LOGBIG').

Is the matrix positive definite (a correlation or covariance matrix)? If so, you can compute the log-det as follows:

Let G = root(A) be the Cholesky root of the matrix A.

Then log(det(A)) = log(det(G`*G)) = log(det(G`)*det(G)) = 2*log(det(G))

Since G is triangular, det(G) = prod(diag(G)) and therefore log(det(G))=sum(diag(G)).

So your formula for log(det(A)) is 2*sum(log(diag(root(A)))). In SAS/IML you would use the VECDIAG function to get the diagonal elements of a matrix.

Message was edited by: Rick Wicklin

Re: EM-algorithm

renc — Wed, 24 Oct 2012 15:05:54 GMT

Thank you, Rick. The idea is very good and helpful. The formula might be log(det(A))= 2*log(det(root(A))). If det(root(A)) is still overflow, we may need to use all eigenvalues of A to calculate the log(det(A)).

Re: EM-algorithm

Rick_SAS — Wed, 24 Oct 2012 15:14:15 GMT

You're free to try it, but computing det(root(A)) is only a minor improvement over det(A). The whole point is to get rid of det(LargeMatrix).

Re: EM-algorithm

Rick_SAS — Wed, 24 Oct 2012 16:55:42 GMT

Rats, I omitted a 'log' in the final formula, which is probably what confused you. I will edit my original post to correct it.