[Genabel-commits] r1658 - pkg/ProbABEL/src

[noreply at r-forge.r-project.org]

Author: lckarssen
Date: 2014-03-21 23:57:20 +0100 (Fri, 21 Mar 2014)
New Revision: 1658

Modified:
   pkg/ProbABEL/src/reg1.cpp
Log:
Speed-ups in ProbABEL's logistic regression.
Profiling showed lots of (expensive) calls to exp() and log(). I got rid of ~ 1/3 of the exp() calls by saving the result in a variable and reusing it in the calculation of exp(mu) / ( 1+exp(mu) ).
The number of log() calls was reduced even more by using the fact that sum_i( log(x_i) ) = log( prod_i(x_i) )

In total this gives roughly a speed up of 30% -- 40% when reading txt dosage files and roughly 20% -- 30% when using filevector files (measured for dosage data).


Modified: pkg/ProbABEL/src/reg1.cpp
===================================================================
--- pkg/ProbABEL/src/reg1.cpp	2014-03-21 21:05:44 UTC (rev 1657)
+++ pkg/ProbABEL/src/reg1.cpp	2014-03-21 22:57:20 UTC (rev 1658)
@@ -719,7 +719,8 @@
             double emu = eMu.get(i, 0);
             double value = emu;
             double zval;
-            value = exp(value) / (1. + exp(value));
+            double expval = exp(value);
+            value = expval / (1. + expval);
             residuals[i] = (rdata.Y).get(i, 0) - value;
             eMu.put(value, i, 0);
             W.put(value * (1. - value), i, 0);
@@ -778,11 +779,24 @@
             beta.print();
         }
         // std::cout << "beta:\n"; beta.print();
-        // compute likelihood
+
+        // Compute the likelihood. The following commented code gives
+        // the 'easy to understand' algorithm. The code that's
+        // actually used is mathematically equivalent (remember:
+        // log(a*b) = log(a)+log(b)), but faster because log() is
+        // relatively expensive.
+        //    for (int i = 0; i < eMu.nrow; i++) {
+        //          loglik += rdata.Y[i] * eMu_us[i] - log(1. +
+        //                    exp(eMu_us[i]));
+        //    }
         prevlik = loglik;
         loglik = 0.;
-        for (int i = 0; i < eMu.nrow; i++)
-            loglik += rdata.Y[i] * eMu_us[i] - log(1. + exp(eMu_us[i]));
+        double logterm = 1;
+        for (int i = 0; i < eMu.nrow; i++) {
+            loglik  += rdata.Y[i] * eMu_us[i];
+            logterm *= 1. + exp(eMu_us[i]);
+        }
+        loglik += - log(logterm);
 
         delta = fabs(1. - (prevlik / loglik));
         niter++;