[Rd] R-2.15.2 changes in computation speed. Numerical precision?

Paul Johnson pauljohn32 at gmail.com
Wed Dec 12 19:14:54 CET 2012

Speaking of optimization and speeding up R calculations...

I mentioned last week I want to speed up calculation of generalized
inverses. On Debian Wheezy with R-2.15.2, I see a huge speedup using a
souped up generalized inverse algorithm published by

V. N. Katsikis, D. Pappas, Fast computing of theMoore-Penrose inverse
matrix, Electronic Journal of Linear Algebra,
17(2008), 637-650.

I was so delighted to see the computation time drop on my Debian
system that I boasted to the WIndows users and gave them a test case.
They answered back "there's no benefits, plus Windows is faster than

That sent me off on a bit of a goose chase, but I think I'm beginning
to understand the situation.  I believe R-2.15.2 introduced a tighter
requirement for precision, thus triggering longer-lasting calculations
in many example scripts. Better algorithms can avoid some of that
slowdown, as you see in this test case.

Here is the test code you can run to see:


It downloads a data file from that same directory and then runs some
multiple imputations with the Amelia package.

Here's the output from my computer


That includes the profile of the calculations that depend on the
ordinary generalized inverse algorithm based on svd and the new one.

See? The KP algorithm is faster.  And just as accurate as
Amelia:::mpinv or MASS::ginv (for details on that, please review my
notes in http://pj.freefaculty.org/scraps/profile/qrginv.R).

So I asked WIndows users for more detailed feedback, including
sessionInfo(), and I noticed that my proposed algorithm is not faster
on Windows--WITH OLD R!

Here's the script output with R-2.15.0, shows no speedup from the
KPginv algorithm


On the same machine, I updated to R-2.15.2, and we see the same
speedup from the KPginv algorithm


After that, I realized it is an R version change, not an OS
difference, I was a bit relieved.

What causes the difference in this case?  In the Amelia code, they try
to avoid doing the generalized inverse by using the ordinary solve(),
and if that fails, then they do the generalized inverse. In R 2.15.0,
the near singularity of the matrix is ignored, but not in R 2.15.2.
The ordinary solve is failing almost all the time, thus triggering the
use of the svd based generalized inverse.  Which is slower.

The Katsikis and Pappas 2008 algorithm is the fastest one I've found
after translating from Matlab to R.  It is not so universally
applicable as svd based methods, it will fail if there are linearly
dependent columns. However, it does tolerate columns of all zeros,
which seems to be the problem case in the particular application I am

I tried very hard to get the newer algorithm described here to go as
fast, but it is way way slower, at least in the implementations I
##  KPP
## Vasilios N. Katsikis, Dimitrios Pappas, Athanassios Petralias.  "An
improved method for
## the computation of the Moore Penrose inverse matrix," Applied
## Mathematics and Computation, 2011

The notes on that are in the qrginv.R file linked above.

The fact that I can't make that newer KPP algorithm go faster,
although the authors show it can go faster in Matlab, leads me to a
bunch of other questions and possibly the need to implement all of
this in C with LAPACK or EIGEN or something like that, but at this
point, I've got to return to my normal job.  If somebody is good at
R's .Call interface and can make a pure C implementation of KPP.

I think the key thing is that with R-2.15.2, there is an svd-related
bottleneck in the multiple imputation algorithms in Amelia. The
replacement version of the function Amelia:::mpinv does reclaim a 30%
time saving, while generating imputations that are identical, so far
as i can tell.


Paul E. Johnson
Professor, Political Science      Assoc. Director
1541 Lilac Lane, Room 504      Center for Research Methods
University of Kansas                 University of Kansas
http://pj.freefaculty.org               http://quant.ku.edu

More information about the R-devel mailing list