[Rcpp-devel] Using sugar expressions for evaluating deviance residuals

Douglas Bates bates at stat.wisc.edu
Sun Aug 15 17:34:47 CEST 2010


On Sun, Aug 15, 2010 at 3:39 AM, Romain Francois
<romain at r-enthusiasts.com> wrote:
> Hi,
>
> Could you provide some minimal example using inline, etc ...

I enclose a comparison of several different versions that I wrote.
Interestingly the old code, called by fx3, does reasonably well,
although I'm not sure if this is a fair comparison as the levels of
optimization by the compiler may not be the same for the code
generated by inline and for the code that is part of R.

As we have seen in the past, extracting the pointers and calling
scalar functions on individual elements does well.  The direct use of
the sugar functions is much easier to understand in the code (it is
essentially the same as writing R expressions) except that ifelse() is
expensive in R but not as a sugar expression.

> Here are a few hints in the meantime. In the next version of Rcpp, I've
> added a few things that help generating sugar functions. This is how for
> example the sugar version of choose is currently implemented :
>
> SUGAR_BLOCK_2(choose      , ::Rf_choose     )
>
> The SUGAR_BLOCK_2 macro (perhaps I should prefix it with RCPP_) lives in
> Rcpp/sugar/SugarBlock.h
>
> You can promote y_log_y to a sugar function like this :
>
> inc <- '
> double y_log_y(double y, double mu){
>    return (y) ? (y * log(y/mu)) : 0;
> }
> SUGAR_BLOCK_2( y_log_y , ::y_log_y )
> '
>
> fx <- cxxfunction( signature( y = "numeric", mu = "numeric" ), '
>        NumericVector res = y_log_y(
>                NumericVector(y),
>                NumericVector(mu)
>        ) ;
>        return res ;
> ', plugin = "Rcpp", includes = inc )
> fx( runif(11) , seq(0, 1, .1 ) )
>
>
> This gives you 3 Rcpp::y_log_y functions :
> - one for when y is a double and mu is a NumericVector (or any numeric sugar
> expression)
> - one for when y is a NumericVector and mu is a double
> - one for both are NumericVector.
>
> However, at the moment it does not take care about recycling so it is up to
> the user to make sure that y and mu have the same length. I'm currently
> thinking about how to deal with recycling.
>
>
>
> rep is a possibility here, but consider this hint in the TODO file:
>
>    o   not sure rep should be lazy, i.e. rep( x, 4 ) fetches x[i] 4 times,
>        maybe we should use LazyVector like in outer to somehow cache the
>        result when it is judged expensive to calculate
>
>
>
> One other thing about sugar is that it has to do many checks for missing
> values so if (with the one defined above) you did call this expression:
>
> 2 * wt * (y_log_y(y, mu) + y_log_y(1.-y, 1.-mu))
>
> (and again, currently binary operators +,*,.. don't take care of recycling)
>
> you would have many checks for missing values. That is fine if you may have
> missing values, because it will propagate them correctly, but it can slow
> things down because they are tested for at every step.
>
> If you are sure that y and mu don't contain missing values, then perhaps one
> thing we can do is embed that information in the data so that sugar does not
> have to check for missing values because it just assumes there are not any.
> Most of the code in sugar contains version that ignores missing values,
> controlled by a template parameter. For example seq_len creates a sugar
> expression where we know for sure that it does not contain missing values.
>
>
> One way perhaps to short circuit this is to first write the code with
> double's and then promote it to sugar:
>
> double resid( double y, double mu, double w){
>        return 2 * wt * (y_log_y(y, mu) + y_log_y(1.-y, 1.-mu) ;
> }
>
> But then you need someone to write SUGAR_BLOCK_3 or write it manually.
>
>
> Another idea would be to have something like NumericVector::import_transform
> but that would take 3 vector parameters instead of 1.
>
>
> Sorry if this email is a bit of a mess, I sort of wrote the ideas as they
> came.
>
> Romain
>
> Le 14/08/10 20:28, Douglas Bates a écrit :
>>
>> In profiling code for generalized linear mixed models for binary
>> responses I found that a substantial portion of the execution time was
>> being spend in evaluating the R functions for inverse links and for
>> the deviance residuals.  As I result I wrote C code (in
>> $RSRC/src/library/stats/src/family.c) or some of those.
>>
>> The way that some research is going I expect that I will soon be in
>> the position where I need to evaluate such expressions even more
>> frequently so it is worthwhile to me to tune it up.
>>
>> In general there will be three NumericVector objects -- y (observed
>> responses), wt (weights) and eta (linear predictors) -- plus an
>> IntegerVector ind.  y, wt and ind will all have length n.  eta's
>> length can be a multiple of n.
>>
>> The index vector, ind, is a factor with k<  n levels so all of its
>> elements are between 1 and k.
>>
>> The objective is to apply the inverse link function to eta, producing
>> the predicted mean response, mu, evaluate the deviance residuals for
>> y, mu and wt and sum the deviance residuals according to the indices
>> in ind.
>>
>> The simplest inverse link is for the logit link
>>
>>  NumericVector mu = 1./(1. + exp(-eta));
>>
>> For the deviance residuals I defined an helper function
>>
>> static R_INLINE
>> double y_log_y(double y, double mu)
>> {
>>     return (y) ? (y * log(y/mu)) : 0;
>> }
>>
>> which has the appropriate limiting behavior when y is zero.  Using
>> that, the deviance residuals can be evaluated as
>>
>> 2 * wt * (y_log_y(y, mu) + y_log_y(1.-y, 1.-mu))
>>
>> That last expression could have different lengths of y and mu but I
>> could use rep to extend y, wt and ind to the desired length.
>>
>> The conditional evaluation in y_log_y is not something that can be
>> ignored because, in many cases y consists only of 0's and 1's so
>> either y_log_y(y, mu) or y_log_y(1.-y, 1.-mu) will fail the condition
>> in the ? expression.
>>
>> Are there suggestions on how best to structure the calculation to make
>> it blazingly fast?
>
>
> --
> Romain Francois
> Professional R Enthusiast
> +33(0) 6 28 91 30 30
> http://romainfrancois.blog.free.fr
> |- http://bit.ly/bzoWrs : Rcpp svn revision 2000
> |- http://bit.ly/b8VNE2 : Rcpp at LondonR, oct 5th
> `- http://bit.ly/aAyra4 : highlight 0.2-2
>
>
-------------- next part --------------
library(Rcpp)
library(inline)
library(rbenchmark)
incl0 <- '
inline double y_log_y(double y, double mu) {
    return y ? y*log(y/mu) : 0;
}
'
code0 <- '
    NumericVector y(yy), eta(ee), wt(ww);
    int n = y.size();
    NumericVector mu = 1./(1. + exp(-eta)), y1 = 1. - y;
    NumericVector yly0(n), yly1(n), mu1 = 1. - mu;
    std::transform(y.begin(), y.end(), mu.begin(), yly0.begin(), *y_log_y);
    std::transform(y1.begin(), y1.end(), mu1.begin(), yly1.begin(), *y_log_y);
    return wrap(2. * wt * (yly0 + yly1));
'
code1 <- '
    NumericVector y(yy), eta(ee), wt(ww);
    NumericVector mu = 1./(1. + exp(-eta)), y1 = 1. - y;
    return wrap(2. * wt * (ifelse(y != 0., y*log(y/mu), 0.) +
                           ifelse(y1 != 0. , y1*(log(y1/(1.-mu))), 0.)));
'
incl2 <- '
inline NumericVector devResid(NumericVector const& y, NumericVector const& mu, NumericVector const& wt) {
    NumericVector y1 = 1. - y;
    return 2. * wt * (ifelse(y != 0., y*log(y/mu), 0.) +
                      ifelse(y1 != 0. , y1*(log(y1/(1.-mu))), 0.));
}
'
code2 <- '
    NumericVector y(yy), eta(ee), wt(ww);
    NumericVector mu = 1./(1. + exp(-eta));
    return wrap(devResid(y, mu, wt));
'
incl4 <- '
inline double y_log_y(double y, double mu) {
    return y ? y*log(y/mu) : 0.;
}
inline double ymu(double y, double mu) {
    return y_log_y(y, mu) + y_log_y(1. - y, 1. - mu);
}
inline double linkinv(double eta) {
    return 1./(1. + exp(-eta));
}
inline double devResid(double y, double mu, double wt) {
    return 2. * wt * ymu(y, mu);
}
'
code4 <- '
    NumericVector y(yy), eta(ee), wt(ww);
    int n = y.size();
    NumericVector ans(n);
    double *ap = ans.begin(), *ep = eta.begin(),
        *wp = wt.begin(), *yp = y.begin();
    for (int i = 0; i < n; i++)
        ap[i] = devResid(yp[i], linkinv(ep[i]), wp[i]);
    return wrap(ans);
'
fx0 <- cxxfunction(signature(yy="numeric", ee="numeric", ww="numeric"),
                   code0, plugin="Rcpp", includes=incl0)
fx1 <- cxxfunction(signature(yy="numeric", ee="numeric", ww="numeric"),
                   code1, plugin="Rcpp")
fx2 <- cxxfunction(signature(yy="numeric", ee="numeric", ww="numeric"),
                   code2, plugin="Rcpp", includes=incl2)
bb <- binomial()
fx3 <- function(yy,ee,ww) bb$dev.resids(yy, bb$linkinv(ee), ww)
fx4 <- cxxfunction(signature(yy="numeric", ee="numeric", ww="numeric"),
                   code4, plugin="Rcpp", includes=incl4)

N <- 38798
set.seed(123454321)
ee <- rnorm(N)
yy <- rbinom(N, 1, 19354/N)
ww <- rep(1, N)
all.equal(fx0(yy,ee,ww), fx1(yy,ee,ww))
all.equal(fx0(yy,ee,ww), fx2(yy,ee,ww))
all.equal(fx0(yy,ee,ww), fx3(yy,ee,ww))
all.equal(fx0(yy,ee,ww), fx4(yy,ee,ww))
benchmark(fx0(yy,ee,ww), fx1(yy,ee,ww),
          fx2(yy,ee,ww), fx3(yy,ee,ww), fx4(yy,ee,ww),
          columns=c("test", "replications", "elapsed", "user.self"))
-------------- next part --------------

R version 2.11.1 (2010-05-31)
Copyright (C) 2010 The R Foundation for Statistical Computing
ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(Rcpp)
> library(inline)
> library(rbenchmark)
> incl0 <- '
+ inline double y_log_y(double y, double mu) {
+     return y ? y*log(y/mu) : 0;
+ }
+ '
> code0 <- '
+     NumericVector y(yy), eta(ee), wt(ww);
+     int n = y.size();
+     NumericVector mu = 1./(1. + exp(-eta)), y1 = 1. - y;
+     NumericVector yly0(n), yly1(n), mu1 = 1. - mu;
+     std::transform(y.begin(), y.end(), mu.begin(), yly0.begin(), *y_log_y);
+     std::transform(y1.begin(), y1.end(), mu1.begin(), yly1.begin(), *y_log_y);
+     return wrap(2. * wt * (yly0 + yly1));
+ '
> code1 <- '
+     NumericVector y(yy), eta(ee), wt(ww);
+     NumericVector mu = 1./(1. + exp(-eta)), y1 = 1. - y;
+     return wrap(2. * wt * (ifelse(y != 0., y*log(y/mu), 0.) +
+                            ifelse(y1 != 0. , y1*(log(y1/(1.-mu))), 0.)));
+ '
> incl2 <- '
+ inline NumericVector devResid(NumericVector const& y, NumericVector const& mu, NumericVector const& wt) {
+     NumericVector y1 = 1. - y;
+     return 2. * wt * (ifelse(y != 0., y*log(y/mu), 0.) +
+                       ifelse(y1 != 0. , y1*(log(y1/(1.-mu))), 0.));
+ }
+ '
> code2 <- '
+     NumericVector y(yy), eta(ee), wt(ww);
+     NumericVector mu = 1./(1. + exp(-eta));
+     return wrap(devResid(y, mu, wt));
+ '
> incl4 <- '
+ inline double y_log_y(double y, double mu) {
+     return y ? y*log(y/mu) : 0.;
+ }
+ inline double ymu(double y, double mu) {
+     return y_log_y(y, mu) + y_log_y(1. - y, 1. - mu);
+ }
+ inline double linkinv(double eta) {
+     return 1./(1. + exp(-eta));
+ }
+ inline double devResid(double y, double mu, double wt) {
+     return 2. * wt * ymu(y, mu);
+ }
+ '
> code4 <- '
+     NumericVector y(yy), eta(ee), wt(ww);
+     int n = y.size();
+     NumericVector ans(n);
+     double *ap = ans.begin(), *ep = eta.begin(),
+         *wp = wt.begin(), *yp = y.begin();
+     for (int i = 0; i < n; i++)
+         ap[i] = devResid(yp[i], linkinv(ep[i]), wp[i]);
+     return wrap(ans);
+ '
> fx0 <- cxxfunction(signature(yy="numeric", ee="numeric", ww="numeric"),
+                    code0, plugin="Rcpp", includes=incl0)
> fx1 <- cxxfunction(signature(yy="numeric", ee="numeric", ww="numeric"),
+                    code1, plugin="Rcpp")
> fx2 <- cxxfunction(signature(yy="numeric", ee="numeric", ww="numeric"),
+                    code2, plugin="Rcpp", includes=incl2)
> bb <- binomial()
> fx3 <- function(yy,ee,ww) bb$dev.resids(yy, bb$linkinv(ee), ww)
> fx4 <- cxxfunction(signature(yy="numeric", ee="numeric", ww="numeric"),
+                    code4, plugin="Rcpp", includes=incl4)
> 
> N <- 38798
> set.seed(123454321)
> ee <- rnorm(N)
> yy <- rbinom(N, 1, 19354/N)
> ww <- rep(1, N)
> all.equal(fx0(yy,ee,ww), fx1(yy,ee,ww))
[1] TRUE
> all.equal(fx0(yy,ee,ww), fx2(yy,ee,ww))
[1] TRUE
> all.equal(fx0(yy,ee,ww), fx3(yy,ee,ww))
[1] TRUE
> all.equal(fx0(yy,ee,ww), fx4(yy,ee,ww))
[1] TRUE
> benchmark(fx0(yy,ee,ww), fx1(yy,ee,ww),
+           fx2(yy,ee,ww), fx3(yy,ee,ww), fx4(yy,ee,ww),
+           columns=c("test", "replications", "elapsed", "user.self"))
             test replications elapsed user.self
1 fx0(yy, ee, ww)          100   1.591      1.56
2 fx1(yy, ee, ww)          100   1.690      1.68
3 fx2(yy, ee, ww)          100   1.631      1.61
4 fx3(yy, ee, ww)          100   0.852      0.83
5 fx4(yy, ee, ww)          100   0.804      0.80
> 
> 
> proc.time()
   user  system elapsed 
 19.590   1.380  21.105 


More information about the Rcpp-devel mailing list