Solution is sensitive to random seed when low-rank approximation is used

[pcarbo]

Here is an example:

library(ashr)
dat <- readRDS("chromium1-ERCC-00054.Rds")
x    <- dat$x
s    <- dat$s
grid <- seq(0,2e-5,length.out = 100)
set.seed(1)
out1 <- ash_pois(x,s,mixcompdist = "halfuniform")
set.seed(2)
out2 <- ash_pois(x,s,mixcompdist = "halfuniform")
print(out1$loglik,digits = 12)
print(out2$loglik,digits = 12)
# -219.299427892
# -140.139152679

The sessionInfo producing this result:

> sessionInfo()
R version 3.6.2 (2019-12-12)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Catalina 10.15.5

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib

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.2-51

loaded via a namespace (and not attached):
 [1] compiler_3.6.2    Matrix_1.2-18     mixsqp_0.3-44     Rcpp_1.0.3
 [5] SQUAREM_2017.10-1 grid_3.6.2        truncnorm_1.0-8   irlba_2.3.3
 [9] invgamma_1.1      lattice_0.20-38

chromium1-ercc-00054.rds.gz

[aksarkar]

The example can be made more concrete by fixing the mode, so that mixsqp is called only once, e.g.,

set.seed(1)
print(ash_pois(x, s, mixcompdist="halfuniform", mode=6e-5, outputlevel="loglik", control=list(verbose=TRUE))$loglik)
set.seed(2)
print(ash_pois(x, s, mixcompdist="halfuniform", mode=6e-5, outputlevel="loglik", control=list(verbose=TRUE))$loglik)