[R-pkg-devel] help with BLAS/LAPACK version debugging

[Iñaki Ucar]

Hi Ben,

I see:

$ FLEXIBLAS=openblas-serial Rscript test.R
Loading required package: Matrix
boundary (singular) fit: see help('isSingular')
Error in solve.default(chol(Hessian), gradient) :
 system is computationally singular: reciprocal condition number = 0
Calls: with ... with.default -> eval -> eval -> solve -> solve -> solve.default

$ FLEXIBLAS=netlib Rscript test.R
Loading required package: Matrix
boundary (singular) fit: see help('isSingular')
Error in h(simpleError(msg, call)) :
 error in evaluating the argument 'a' in selecting a method for
function 'solve': the leading minor of order 2 is not positive
definite
Calls: with ... chol -> chol -> chol.default -> .handleSimpleError -> h

This is for Fedora 35 with R 4.1.3 and FlexiBLAS 3.1.3, which allows
you to easily switch (as shown above) between several implementations
available, including NETLIB 3.10.0, OpenBLAS 0.3.19 (serial, openmp,
threads), BLIS 0.7.0 (serial, openmp, threads) and ATLAS 3.10.3, but
also other libraries you already have if you tell FlexiBLAS where they
are. If you want to play with those, more info at [1], [2], and of
course [3].

Iñaki

[1] https://cran.r-project.org/bin/linux/fedora/#blaslapack-switching
[2] https://docs.fedoraproject.org/en-US/packaging-guidelines/BLAS_LAPACK/#_backend_selection
[3] https://github.com/mpimd-csc/flexiblas/#documentation

On Tue, 12 Apr 2022 at 02:34, Ben Bolker <bbolker using gmail.com> wrote:
>
>    Hi folks,
>
>    Apologies for long post. (Part 1, overview; part 2, list of attempts
> so far; part 3, test code.)
>
> OVERVIEW
>
>    CRAN is reporting check errors with the latest release of lme4
> <https://cran.r-project.org/web/checks/check_results_lme4.html> that
> *may* (best guess) be caused by some linear algebra returning NA values
> rather than throwing an error as expected.
>
>    Kurt Hornik says:
>
>    > The failing check runs are using an external LAPACK 3.10.0. To be
> even more precise, Brian does not get this on Fedora with a system
> reference (netlib) 3.10 LAPACK, whereas I do on Debian with OpenBLAS.  Hth.
>
>    The relevant test code is below.
>
>    Further notes:
> https://github.com/lme4/lme4/blob/master/misc/R_lapackstuff.md ;
> https://github.com/lme4/lme4/issues/677
>
>    I think I probably know how to fix/work around the problem, but I
> would dearly love to have a test system where I could (a) confirm that
> things are failing in the way/for the reason I expect; and (b) test that
> my proposed fix actually works.
>
>    Any ideas (or reports that the test code fails on your system) would
> be greatly appreciated ...
>
>    cheers
>      Ben Bolker
> =====
>
> ATTEMPTS SO FAR (none reproduce the error)
>
> * running on r-hub debian-gcc-devel, both remotely and in a local Docker
> container. AFAICT it is running
>
> BLAS:   /opt/R-devel/lib/R/lib/libRblas.so
> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
>
>    I was unsuccessful in using Debian's update-defaults command to mess
> around with the LAPACK version.
>
>   * running on the `rocker/verse` image from the Rocker project.
>
> The default is openblas-pthread/liblapack.so.3; I don't get any errors
> with this. I can use update-defaults to switch to liblapack.so.3.9.0
> (not 3.10.0, this is using Debian stable I think), also no errors that way.
>
> * running on my own machine (PopOS!, r-devel)
>
> This is using openblas-pthread/libopenblasp-r0.3.13.so, and I can't for
> the life of me figure out how to switch BLAS/LAPACK versions (have tried
> reinstalling both with (--with-lapack --with-blas) and without, using
> update-alternatives doesn't seem to make any difference?  Have read
> various online docs about how this is supposed to work ...
>
>
> =======
>
> TEST CODE
>
> ## 1. this succeeds for me, fails for CRAN
>
> set.seed(101)
> d <- data.frame(z=rnorm(200),
>                  f=factor(sample(1:10,200, replace=TRUE)))
> library(lme4)
> library(testthat)
> fm1 <- lmer(z~ as.numeric(f) + 1|f, d)
> fm1 using optinfo$derivs$Hessian[2,2] <- NA
> expect_warning(lme4:::checkConv(fm1 using optinfo\$derivs,
>                                  coefs=c(1,1),
>                                  ctrl=lmerControl()$checkConv,lbound=0),
>                 "Problem with Hessian check")
>
> ## 2. this throws 'system is computationally singular' for me, returns a
> vector containing NAs for CRAN (I think)
>
> dd <- list(gradient = c(0.00132136676711525, 0.00898197413334856, 0
> ), Hessian = structure(c(195.253128051758, 962.483270645142,
> 0, 962.483270645142, NA, 0, 0, 0, 6542.44775390625), dim = c(3L,
> 3L)))
> with(dd, solve(chol(Hessian),gradient))
>
> ## 3. report BLAS/LAPACK info
>
> s <- sessionInfo()
> cat(s$BLAS, s$LAPACK, sep = "\n")
>
> ______________________________________________
> R-package-devel using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-package-devel



-- 
Iñaki Úcar