Hi all,
On Tue, Jan 1, 2013 at 6:01 PM, Thomas Lumley <tlum...@uw.edu> wrote: > Now, if most of the matrices you are trying to invert are actually > invertible (as I would hope), it may be quicker to use the Cholesky > approach will a fallback to the SVD for semidefinite matrices. That is, > something like > > tryCatch( > chol2inv(chol(xx)), > error=function(e) ginv(xx) > ) > > Just to be clear, the Amelia package takes this approach, since our EM algorithm can wander off into questionable parts of the parameter space (with non-invertible matrices). The issue that Paul pointed out was an issue my implementation of the Cholesky approach which broke the try() function and so used SVD for every inverse, which of course is slow(er). This has been fixed in the most recent version (1.6.4) and we are finishing up a version that uses Rcpp to speed things up quite a bit. Cheers, matt. ~~~~~~~~~~~ Matthew Blackwell Assistant Professor of Political Science University of Rochester url: http://www.mattblackwell.org > Most of the time you will get the Cholesky branch, which is much faster > (about five-fold for 10x10 matrices on my system). On my system and using > a 10x10 matrix the overhead in the tryCatch() is much smaller than the time > taken by either set of linear algebra, so there is a net gain as long as > even a reasonably minority of the matrices are actually invertible. > > You can probably do slightly better by replacing the chol2inv() with > backsolve(): solving just the systems of linear equations you need is > usually preferable to constructing a matrix inverse. > > > Note that this approach will give wrong answers without warning if the > matrices are not symmetric. > > -thomas > > -- > Thomas Lumley > Professor of Biostatistics > University of Auckland > > [[alternative HTML version deleted]] > > ______________________________________________ > R-devel@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-devel > [[alternative HTML version deleted]] ______________________________________________ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
