Basic operations in R giving different results on Windows and Linux – Debian

I ran your code on my local machine with Fedora 37 and R 4.2.2. As the other commenters, I also got the result you got on Windows. Then I pulled the rocker image for R version 3.6.3:

docker run -ti rocker/r-ver:3.6.3 R

This image is also Debian-based. Here I was able to recreate the result you got on your system.

Then I moved to the rocker versioned release 4.0.0:

docker run -ti rocker/r-ver:4.0.0 R

Here the result was the same as you got on Windows and everyone else got on their machine. It must be noted that with R 4.0.0 the rocker project moved from Debian-based images to Ubuntu LTS.

Fedora comes with the ability to easily switch the BLAS/LAPACK backend via the flexiblas package. Thanks to that, I was able to test your code with the eight different backends available on my system. As you can see below, they do yield different results. In particular, the ATLAS backend comes somewhat close to the result you got. In contrast, OPENBLAS-OPENMP (the default on Fedora), other OPENBLAS variants, and NETLIB all produce the same result as you received on Windows. A third family BLIS produces yet another set of possible results.

Is one of the results better than the others? Yes! optim() looks for a result that minimizes the supplied function. In its returned list, it reports not just the minimizing parameters, but also the value for them. I’ve included that in the table below. So the ATLAS backend wins the prize here. It must be said that optim() does NOT minimize analytically. So it always gives approximate results. That is why initial values and the method matter for what results we get. And apparently, with your function the backend also matters. And if you look at the parameters you got on Buster the function goes to 0.002800321. So it is actually a better result than what we all get on our more modern systems, except for the result I got with ATLAS. That also happens to be much slower than the other backends. So it seems, the newer backends might have traded speed for accuracy.

If your aim is consistency across platforms, you can upgrade your system to Debian 11 Bullseye, since that appears to have a backend producing the same results as other modern platforms, as the answer by @jay.sf indicates. You could also investigate if you can find the same BLAS backend version used on Buster for Windows.

Furthermore, you can try to change to another blas library on your current system. Here is a guide how to do that. Though it is for Ubuntu, as both use apt, it should work for your system as well. (Edit: I tried that in a VM for Buster. None of the available BLAS backends produced the same result as on the more modern systems)

Finally, if you feel you must have a newer BLAS library on your older system, then you could try to backport it yourself. I have no experience with this. I don’t know how advisable it is or how likely to succeed. I am just mentioning it for completeness.

library(flexiblas)
library(tidyverse)

test_fun <- function(i) {
  flexiblas_switch(i)
  data <- structure(list(x = c(0.1, 0.1, 0.1, 5, 5, 5, 10, 10, 10, 20, 20, 20),
                         y = c(0.013624804, 0.014023006, 0.013169554, 0.70540352,
                               0.68711807, 0.69233506, 1.4235181, 1.348244, 1.4141854, 2.779813,
                               2.7567347, 2.7436437)), class = c("data.frame"), row.names = c(NA, 12L))
  val <- c(43.3065849160736, 0.00134925463859564, 1.03218302435548, 270.328323775978)
  theta <- 1.60812569803848
  init <- c(b0 = 2.76836653333333, b1 = 0.0134350095, b2 = 2.15105945932773,
            b3 = 6.85922519794374)
  f <- function(X, b0, b1, b2, b3) {
    b0 + (b1 - b0) / (1 + exp(b2*(log(X) - log(b3))))
  }
  W <- 1 / f(data$x, val[1], val[2], val[3], val[4])^theta
  SSw <- function(Y, X, b0, b1, b2, b3, w) {
    sum(w * (Y - f(X, b0, b1, b2, b3))^2)
  }
  SSw.vec <- function(par) SSw(data$y, data$x, par[1], par[2], par[3], par[4], W)
  mod <- optim(init, SSw.vec, method = "L-BFGS-B", lower = c(-Inf,-Inf,-Inf,0))
  return(c(mod$par, value = mod$value))
}

flexiblas_list() |>
  setdiff("__FALLBACK__") |>
  tibble(backend = _) |>
  mutate(
    idx = flexiblas_load_backend(backend),
    res = map(idx, test_fun)
  ) |>
  unnest_wider(res)
#> # A tibble: 8 × 7
#>   backend            idx    b0      b1    b2    b3   value
#>   <chr>            <int> <dbl>   <dbl> <dbl> <dbl>   <dbl>
#> 1 NETLIB               2  31.0 0.00183  1.05  184. 0.00282
#> 2 OPENBLAS-OPENMP      1  31.0 0.00183  1.05  184. 0.00282
#> 3 ATLAS                3  34.7 0.00168  1.04  209. 0.00280
#> 4 BLIS-SERIAL          4  27.1 0.00225  1.06  158. 0.00285
#> 5 BLIS-OPENMP          5  27.1 0.00225  1.06  158. 0.00285
#> 6 BLIS-THREADS         6  27.1 0.00225  1.06  158. 0.00285
#> 7 OPENBLAS-SERIAL      7  31.0 0.00183  1.05  184. 0.00282
#> 8 OPENBLAS-THREADS     8  31.0 0.00183  1.05  184. 0.00282