solve(): system is singular

[william-denault]

Hello,

I am using mixsqp for a project. However, I notice a problem that happens (especially on a cluster) where mixsqp struggles to find a simple solution. Please find a minimal example that runs well on a laptop but creates problems on clusters:

library(ashr)
library(mixsqp)

bethat <- c(0.123762324,
            0.097227906,
            -0.029804568,
            0.020958705,
            -0.002632602,
            -0.007505082,
            0.097227906,
            -0.034825676,
            -0.034825676,
            -0.050515880)

sdhat <- c(0.04338368,
           0.04052425,
           0.07405885,
           0.02930913,
           0.03372245,
           0.03956477,
           0.04052425,
           0.02597307,
           0.02597307,
           0.10213671)

sd_list <- c(0.000000000,
             0.002561258,
             0.003622166,
             0.005122516,
             0.007244332,
             0.010245033,
             0.014488664,
             0.020490066,
             0.028977329,
             0.040980131,
             0.057954657,
             0.081960262,
             0.115909314,
             0.163920524,
             0.231818629)

fitash <- ash(bethat,sdhat,mixcompdist = "normal")

sdmat <- sqrt(outer(c(sdhat^2),sd_list^2,"+"))
L <- (dnorm(outer(c(bethat),rep(0,length(sd_list)),FUN = "-")/sdmat,
            log = TRUE) - log(sdmat))

# Dealing in case of due to small sd due to small sample size.
L <- apply(L,2,function (x) {
    x[which(is.na(x))] <- median(x,na.rm = TRUE)
    return(x)
  })

out <- mixsqp(L,log = TRUE)

The error I got :

`warning: solve(): system is singular (rcond: 6.21885e-18); attempting approx solution

warning: solve(): system is singular (rcond: 1.49366e-16); attempting approx solution

warning: solve(): system is singular (rcond: 1.55578e-16); attempting approx solution

warning: solve(): system is singular (rcond: 7.29899e-18); attempting approx solution

warning: solve(): system is singular (rcond: 1.79796e-16); attempting approx solution

warning: solve(): system is singular (rcond: 2.00388e-16); attempting approx solution`

current session parameter on cluster

` R version 4.1.0 (2021-05-18) Platform: x86_64-pc-linux-gnu (64-bit) Running under: CentOS Linux 7 (Core)

Matrix products: default BLAS: /software/R-4.1.0-no-openblas-el7-x86_64/lib64/R/lib/libRblas.so LAPACK: /software/R-4.1.0-no-openblas-el7-x86_64/lib64/R/lib/libRlapack.so

locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=C LC_COLLATE=C LC_MONETARY=C LC_MESSAGES=C LC_PAPER=C LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=C LC_IDENTIFICATION=C

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

other attached packages: [1] susiF.alpha_0.1.161 mixsqp_0.3-43 ashr_2.2-54 wavethresh_4.7.0 MASS_7.3-54

loaded via a namespace (and not attached): [1] compiler_4.1.0 Matrix_1.3-3 tools_4.1.0 Rcpp_1.0.9 SQUAREM_2021.1 grid_4.1.0 truncnorm_1.0-8 irlba_2.3.5 invgamma_1.1 lattice_0.20-44`

[pcarbo]

@william-denault What version of RcppArmadillo do you have?

[william-denault]

@pcarbo I have the latest version (version 0.11.2.4.0)

[pcarbo]

I can confirm that I get the solution just fine on my laptop as well (an Intel MacBook Pro).

[pcarbo]

@william-denault I notice you are not using the latest versions of mixsqp and ashr. Can you please try again with the latest versions?

[william-denault]

@pcarbo I have updated the packages and the problem persists with quite a substantial difference in terms of estimated pi.

[pcarbo]

Any particular reason you are using R without openblas? What happens when you try the R with openblas?

[william-denault]

The problem remains, unfortunately, yet the solutions seems quite similar to what I have on my laptop but I still have warnings about matrix conditionning

[pcarbo]

It is reassuring that you get a similar solution. You may be able to avoid these warnings by adjusting some of the control parameters, e.g., tol.svd, eps, or you can simply ignore the warning and move on.