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