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

[Paul Johnson]

On Thu, Dec 13, 2012 at 3:33 AM, Uwe Ligges
<ligges at statistik.tu-dortmund.de> wrote:
> Long message, but as far as I can see, this is not about base R but the
> contributed package Amelia: Please discuss possible improvements with its
> maintainer.
>
Thanks for answering, but I'm really surprised by your answer. The
package is a constant, but R got much slower  between R-2.15.1 and
R-2.15.2. The profile of functions called radically changed, svd gets
called much more because solve() fails much more often.

No change in R could account for it?  I'm not saying R is wrong, it
may be more accurate and better. After chasing the slowdown for a
week, I need some comfort. Does the LINPACK -> LAPACK change play a
role. The change I'm looking for is something that would substantially
tune up mathematical precision so that matrices that did not seem
singular before are now, in the eyes of functions like chol() and
solve().  Whereas in R-2.15.1, a matrix can be inverted by solve(),
for example, now R-2.15.2 rejects the matrix as singular.

I will take the problem up with applications, of course. But you see
how package writers might think its ridiculous.   They argue, "I had a
perfectly fine calculation against R-2.15.0 and R-2.15.1, and now with
R-2.15.2 it takes three times as long?  And you want me to revise my
package?"

Would you be persuaded there's an R base question if I showed you a
particular matrix that can be decomposed or solved in R-2.15.1 but
cannot be in R-2.15.2? I should have thought of that before, I suppose
:)

pj
> Best,
> Uwe Ligges
>
>
>
> On 12.12.2012 19:14, Paul Johnson wrote:
>>
>> 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
>> Linux".
>>
>> 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:
>>
>> http://pj.freefaculty.org/scraps/profile/prof-puzzle-1.R
>>
>> 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
>>
>> http://pj.freefaculty.org/scraps/profile/prof-puzzle-1.Rout
>>
>> 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
>>
>> http://pj.freefaculty.org/scraps/profile/prof-puzzle-1-Windows.Rout
>>
>> On the same machine, I updated to R-2.15.2, and we see the same
>> speedup from the KPginv algorithm
>>
>>
>> http://pj.freefaculty.org/scraps/profile/prof-puzzle-1-CRMDA02-WinR2.15.2.Rout
>>
>> 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
>> testing.
>>
>> 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
>> tried:
>> ##  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.
>>
>> pj
>>
>>
>>
>



-- 
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