Scaling of mixIP with n and k

Overview

Here we assess how the convex optimization in mixIP scales with n and k.

mixIP: Number of samples

#define a function that calls the basic function to do mixIP in ash
mixIP = function(lik,rtol=1e-6){
  k = ncol(lik)
  prior=rep(1,k)
  prior[1]=10
  n = nrow(lik)
  init = ashr:::initpi(k,n,1)
  ashr::mixIP(lik,prior,init,ashr:::set_control_mixIP(list(rtol=rtol)))
}

Idea is to get a quick idea of how mixIP scales with number of samples:

  timer_n = function(seed){
    set.seed(seed)
    z = rnorm(100000,0,2)
    #gm = exp(seq(0.01,1,length=10)) # grid of multiplier
    ll = ashr::ash(z,1,mixcompdist= "normal",outputlevel=4)
    lik = ll$fit_details$matrix_lik
    
    time = list()
    time[[1]] = system.time(mixIP(lik))
    time[[2]] = system.time(mixIP(lik[1:10000,]))
    time[[3]] = system.time(mixIP(lik[1:1000,]))
    return(time)
  }
  timer_n(1)

[[1]]
   user  system elapsed 
 14.805   2.397   7.969 

[[2]]
   user  system elapsed 
  1.242   0.243   0.579 

[[3]]
   user  system elapsed 
  0.099   0.019   0.037 

  timer_n(2)

[[1]]
   user  system elapsed 
 17.697   2.757   7.354 

[[2]]
   user  system elapsed 
  1.364   0.256   0.593 

[[3]]
   user  system elapsed 
  0.096   0.021   0.039 

quick conclusion: time is more than linear. We might benefit from, say, averaging results over subsets.

gridsizes

 timer_k = function(seed){
    set.seed(seed)
    z = rnorm(10000,0,2)
    #gm = exp(seq(0.01,1,length=10)) # grid of multiplier
    ll = ashr::ash(z,1,mixcompdist= "normal", gridmult = 1.05, outputlevel=4)
    lik = ll$fit_details$matrix_lik
     time = list()
     k = ncol(lik)
    time[[1]] = system.time(mixIP(lik))
    time[[2]] = system.time(mixIP(lik[,c(1,sample(2:k,round(k/2),replace=F))]))
    time[[3]] = system.time(mixIP(lik[,c(1,sample(2:k,round(k/4),replace=F))]))
    time[[4]] = system.time(mixIP(lik[,c(1,sample(2:k,round(k/8),replace=F))]))
    return(time)
  }
  timer_k(1)

[[1]]
   user  system elapsed 
  4.860   0.485   1.865 

[[2]]
   user  system elapsed 
  2.305   0.280   1.005 

[[3]]
   user  system elapsed 
  1.766   0.353   0.726 

[[4]]
   user  system elapsed 
  1.032   0.209   0.436 

  timer_k(2)

[[1]]
   user  system elapsed 
  4.709   0.550   1.938 

[[2]]
   user  system elapsed 
  3.034   0.396   1.155 

[[3]]
   user  system elapsed 
  1.598   0.345   0.665 

[[4]]
   user  system elapsed 
  0.930   0.200   0.397 

Rough conclusion: halving k almost halves time (but not quite), so we might do well to use dense grid on small n to get plausible grid values, and then narrow down for larger n.

Session information

sessionInfo()

R version 3.3.1 (2016-06-21)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.11.5 (El Capitan)

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] ashr_2.1

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.8       knitr_1.15.1      magrittr_1.5     
 [4] REBayes_0.73      MASS_7.3-45       doParallel_1.0.10
 [7] pscl_1.4.9        SQUAREM_2016.8-2  lattice_0.20-34  
[10] foreach_1.4.3     stringr_1.1.0     tools_3.3.1      
[13] parallel_3.3.1    grid_3.3.1        htmltools_0.3.5  
[16] iterators_1.0.8   assertthat_0.1    yaml_2.1.14      
[19] rprojroot_1.1     digest_0.6.10     Matrix_1.2-7.1   
[22] codetools_0.2-15  evaluate_0.10     rmarkdown_1.3    
[25] stringi_1.1.2     Rmosek_7.1.2      backports_1.0.4  
[28] truncnorm_1.0-7