How to compute the inverse of a close to singular matrix in R?

Question

I want to minimize function FlogV (working with a multinormal distribution, Z is data matrix NxC; SIGMA it´s a square matrix CxC of var-covariance of data, R a vector with length C)

FLogV <- function(P){

(here I define parameters, P, within R and SIGMA)

logC <- (C/2)*N*log(2*pi)+(1/2)*N*log(det(SIGMA))

SOMA.t <- 0

for (j in 1:N){
SOMA.t <- SOMA.t+sum(t(Z[j,]-R)%*%solve(SIGMA)%*%(Z[j,]-R))
}

MlogV <- logC + (1/2)*SOMA.t

return(MlogV)

}

minLogV <- optim(P,FLogV)

All this is part of an extend code which was already tested and works well, except in the most important thing: I can´t optimize because I get this error: “Error in solve.default(SIGMA) : system is computationally singular: reciprocal condition number = 3.57726e-55”

If I use ginv() or pseudoinverse() or qr.solve() I get: “Error in svd(X) : infinite or missing values in 'x'”

The thing is: if I take the SIGMA matrix after the error message, I can solve(SIGMA), the eigen values are all positive and the determinant is very small but positive det(SIGMA) [1] 3.384674e-76

eigen(SIGMA)$values
[1] 0.066490265 0.024034173 0.018738777 0.015718562 0.013568884 0.013086845
….
[31] 0.002414433 0.002061556 0.001795105 0.001607811

I already read several papers about change matrices like SIGMA (which are close to singular), did several transformations on data scale and form but I realized that, for a 34x34 matrix like the example, after det(SIGMA) close to e-40, R assumes it like 0 and calculation fails; also I can´t reduce matrix dimensions and can´t input in my function correction algorithms to singular matrices because R can´t evaluate it working with this optimization functions like optim. I really appreciate any suggestion to this problem. Thanks in advance, Maria D.

Neal Fultz Neal Fultz · Accepted Answer · 2013-11-20T01:34:33

It isn't clear from your post whether the failure is coming from det() or solve()

If its just the solve in the quadratic term, you may want to try the two argument version of solve, it can be a bit more stable. solve(X,Y) is the same as solve(X) %*% Y

If you can factor sigma using chol(), you will get a triangular matrix such that LL'=Sigma. The determinant is the product of the diagonals, and you might try this for the quadratic term:

crossprod( backsolve(L, Z[j,]-R))

How to compute the inverse of a close to singular matrix in R?

1 Answers