Statistical programming, matrix languages, and more

EM-algorithm

Reply
Occasional Contributor
Posts: 11

EM-algorithm

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!


SAS Super FREQ
Posts: 3,615

Re: EM-algorithm

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

Occasional Contributor
Posts: 11

Re: EM-algorithm

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

SAS Super FREQ
Posts: 3,615

Re: EM-algorithm

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

SAS Super FREQ
Posts: 3,615

Re: EM-algorithm

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

Ask a Question
Discussion stats
  • 4 replies
  • 432 views
  • 1 like
  • 2 in conversation