In Mouse-Imaging-Centre/RMINC: Statistical Tools for Medical Imaging NetCDF (MINC) Files

knitr::opts_chunk$set(echo = TRUE, cache = TRUE)

RcppArmadillo provides a fast lm function implemented in c++, I suspect it would be an easier starting point for getting residuals than hacking through the gorey guts of minc2_model. I want to figure out how much performance we can expect to lose when converting away

We'll start with a small test from test_mincLm

Test 1

Let's load up the test files

suppressPackageStartupMessages({
  library(RMINC)
  library(dplyr)
  library(Rcpp)
  library(RcppArmadillo)
  library(rbenchmark)
})

getRMINCTestData("./")

gf <- read.csv("rminctestdata/test_data_set.csv")

voxel_left <- mincGetVoxel(gf$jacobians_fixed_2[1:10], 0,0,0)
voxel_right <- mincGetVoxel(gf$jacobians_fixed_2[11:20], 0,0,0)
Sex <- gf$Sex[1:10]
Scale <- gf$scale[1:10]
Coil <- as.factor(gf$coil[1:10])

gf$coil <- as.factor(gf$coil)
gftest <- gf[1:10,]

Now we'll setup the rcppArmadillo function, from https://github.com/RcppCore/RcppArmadillo/blob/master/inst/examples/fastLm.r.

src <- '
Rcpp::List fLmSEXP(SEXP Xs, SEXP ys) {
    Rcpp::NumericMatrix Xr(Xs);
    Rcpp::NumericVector yr(ys);
    int n = Xr.nrow(), k = Xr.ncol();
    arma::mat X(Xr.begin(), n, k, false);
    arma::colvec y(yr.begin(), yr.size(), false);
    int df = n - k;
    // fit model y ~ X, extract residuals
    arma::colvec coef = arma::solve(X, y);
    arma::colvec res  = y - X*coef;
    double s2 = std::inner_product(res.begin(), res.end(),
                                   res.begin(), 0.0)/df;
    // std.errors of coefficients
    arma::colvec sderr = arma::sqrt(s2 *
       arma::diagvec(arma::pinv(arma::trans(X)*X)));
    return Rcpp::List::create(Rcpp::Named("coefficients")=coef,
                              Rcpp::Named("stderr")      =sderr,
                              Rcpp::Named("df")          =df);
}
'
cppFunction(code=src, depends="RcppArmadillo") #creates fLmSEXP

Now we'll benchmark them

{# Hide this code from knitr's output capture (it breaks normal sinking)
  sink("/dev/null")
  pred_mat <- as.matrix(as.numeric(gftest$Sex))
  benches <- benchmark(mincLm = mincLm(jacobians_fixed_2 ~ Sex, data = gftest)
                       , fLmSEXP = mincApplyRCPP(gftest$jacobians_fixed_2
                                                 , function(x) fLmSEXP(pred_mat, as.matrix(x))
                                                 , slab_sizes = c(10,10,1))
                       , lm.fit = mincApplyRCPP(gftest$jacobians_fixed_2
                                                , function(x) lm.fit(pred_mat, as.matrix(x))
                                                , slab_sizes = c(10,10,1))
                       , columns = c("test", "replications", "relative", "elapsed")
                       , order = "relative")
  sink()
}

# Runtime
benches

At time of first run this looks like mincApplyRCPP w/ fLmSEXP is about 20% slower than mincLm, and 20% faster than lm.fit. This is not bad given no effort is made to re-use data structures between calls. I wonder how the two would stack up running on a single voxel.

like <- gftest$jacobians_fixed_2[1]
c(1, rep(0, 999)) %>%
  as.minc %>%
  `likeVolume<-`(like) %>%
  mincWriteVolume("one_el_mask.mnc")

{# Hide this code from knitr's output capture (it breaks normal sinking)
  sink("/dev/null")
  benches <- benchmark(mincLm = mincLm(jacobians_fixed_2 ~ Sex, data = gftest, mask = "one_el_mask.mnc")
                       , fLmSEXP = mincApplyRCPP(gftest$jacobians_fixed_2
                                                 , function(x) fLmSEXP(pred_mat, as.matrix(x))
                                                 , slab_sizes = c(10,10,1), mask = "one_el_mask.mnc")
                       , columns = c("test", "replications", "relative", "elapsed")
                       , order = "relative")
  sink()
}

#Runtime
benches

So it looks like for a single voxel, fLmSEXP wins by a few percent. Nothing to write home about. But from a code clarity standpoint, the ARMA code is much nicer. Although I'm a bit unsure of which of the two functions has more overhead. The rcpp version has an additional sort step, and has much larger output objects (list instead of vector), but the standard mincLm has to deal with the formula.

Although mincLm actually does a bit more, it computes the t-statistics. I think the smart thing to do is to piggy-back on minc_voxel_lm.

In addition as Eddelbuettel points out the fit using dqrls is more robust, implementing a better pivoting strategy. Oh well, worth exploring alternatives.

unlink("one_el_mask.mnc")
unlink("rminctestdata", recursive = TRUE)
unlink("rminctestdata.tar.gz")

Mouse-Imaging-Centre/RMINC documentation built on Nov. 12, 2022, 1:50 p.m.