Wednesday 15 September 2010

Comparing matrix inversions in R - what is wrong with the Cholesky method? -



Comparing matrix inversions in R - what is wrong with the Cholesky method? -

i compared various methods compute inverse of symmetric matrix:

solve (from bundle lapck) solve (but using higher machine precision) qr.solve (said faster) ginv (mass package, implementation of moore-penrose algo) chol2inv (using cholesky decomposition)

the inverse-matrix compared through eigenvalues:

r library(mass) ## create matrix m = replicate(10, runif(n=10)) m[lower.tri(m)] = t(m)[lower.tri(m)] ## inverse matrix inv1 = solve(m) inv2 = solve(m, tol = .machine$double.eps) inv3 = qr.solve(m) inv4 = ginv(m) inv5 = chol2inv(m) ## eigenvalues of inverse em1=eigen(inv1) em2=eigen(inv2) em3=eigen(inv3) em4=eigen(inv4) em5=eigen(inv5) ## plot abs of eigenvalues (may complex) mypch=c( 20, 15, 17, 25, 3 ) plot(abs(em1$values),pch=mypch[1],cex=1.5) points(abs(em2$values),pch=mypch[2], cex=1.5) points(abs(em3$values),pch=mypch[3], cex=1.5) points(abs(em4$values),pch=mypch[4], cex=1.5) points(abs(em5$values),pch=mypch[5], cex=1.5) legend( "topright", c("solve","solve-double","solve-fast","moore-penrose","cholesky"), pch=mypch )

as can see inverse given cholesky method different other.

according post, if matrix symmetric (in our case yes), cholesky method preferred: matrix inversion or cholesky?

but solve() beingness "official-wellspread" r method invert method, may rather misunderstand something...

any tip?

thanks in advance,

you need pass cholesky decomposition chol2inv:

inv5 = chol2inv(chol(m))

if m positive definite (which isn't not-reproducible input) should give same result other methods.

r matrix matrix-inverse

No comments:

Post a Comment