Regarding the inability to reproduce results exactly as output, I wrote a script to fit the insteval model allowing for AppleAccelerate to be turned on or off, in which case OpenBLAS will be used. I will check that in to a branch (feel free to move it to another location if you wish).
I tried fitting with OpenBLAS using 1 and 6 threads. The results were the same except for the benchmark timing (which was nearly the same) and the versioninfo() output.
I used the converged parameter estimates from the OpenBLAS runs, called thetaOB1, as values at which to evaluate the objective.
I also ran on an M1-Mac using AppleAccelerate and got slightly different answers
% diff timings_openblas1.jlout timings_accelerate6.jlout
10,11c10,12
< Threads: 1 default, 0 interactive, 1 GC (on 6 virtual cores)
< Status `/private/var/folders/d9/t94kd5m50rl1r3bchncn62sm0000gn/T/jl_MnaLDj/Project.toml`
---
> Threads: 6 default, 0 interactive, 3 GC (on 6 virtual cores)
> Status `/private/var/folders/d9/t94kd5m50rl1r3bchncn62sm0000gn/T/jl_9SSsj9/Project.toml`
> [13e28ba4] AppleAccelerate v0.4.0
17c18
< LBTConfig([ILP64] libopenblas64_.dylib)
---
> LBTConfig([LP64] Accelerate, [ILP64] Accelerate)
33c34
< Function evaluations: 181
---
> Function evaluations: 173
36,37c37,38
< Final parameter vector: [0.27572694783915613, 0.4352917263395991, 0.04316230740526005, 0.1299749675679518]
< Final objective value: 237648.60164780272
---
> Final parameter vector: [0.27572689108074344, 0.4352913075607651, 0.04316427749556987, 0.12998231312662356]
> Final objective value: 237648.60164785833
41c42
< objective(updateL!(setθ!(m1, thetaOB1))) = 237648.60164780272
---
> objective(updateL!(setθ!(m1, thetaOB1))) = 237648.60164780298
43c44
< Sample(time=0.007044291, allocs=58, bytes=1504)
---
> Sample(time=0.005501458000000001, allocs=58, bytes=1504)
I will run this on an Intel-based laptop with OpenBLAS and with MKL for comparison.
My response to the editors is that we are proposing a method to evaluate the profiled log-likelihood for a linear mixed model, which will define the objective to be optimized. The optimization itself will be performed by other software - in our case implementations of Powell's BOBYQA algorithm in Nlopt or PRIMA. Very small changes in the evaluated objective, due to different processor architectures or choices like the BLAS implementation, can result in different optimization trajectories and different converged estimates. The bottom line on whether the differences are important is the difference in the objective at the estimates. Differences are judged according to an analysis of deviance scale which gives an absolute scale indicating that differences on the order of 0.001 can safely be ignored.
So my immediate tasks are to bring PRIMA into the comparisons and to run everything on an Intel-based laptop with and without MKL. Suggestions on presentation of the results and discussion are, of course, welcome.
Regarding the inability to reproduce results exactly as output, I wrote a script to fit the insteval model allowing for AppleAccelerate to be turned on or off, in which case OpenBLAS will be used. I will check that in to a branch (feel free to move it to another location if you wish).
I tried fitting with OpenBLAS using 1 and 6 threads. The results were the same except for the benchmark timing (which was nearly the same) and the
versioninfo()output.I used the converged parameter estimates from the OpenBLAS runs, called
thetaOB1, as values at which to evaluate the objective.I also ran on an M1-Mac using AppleAccelerate and got slightly different answers
I will run this on an Intel-based laptop with OpenBLAS and with MKL for comparison.
My response to the editors is that we are proposing a method to evaluate the profiled log-likelihood for a linear mixed model, which will define the objective to be optimized. The optimization itself will be performed by other software - in our case implementations of Powell's BOBYQA algorithm in
NloptorPRIMA. Very small changes in the evaluated objective, due to different processor architectures or choices like the BLAS implementation, can result in different optimization trajectories and different converged estimates. The bottom line on whether the differences are important is the difference in the objective at the estimates. Differences are judged according to an analysis of deviance scale which gives an absolute scale indicating that differences on the order of 0.001 can safely be ignored.So my immediate tasks are to bring PRIMA into the comparisons and to run everything on an Intel-based laptop with and without MKL. Suggestions on presentation of the results and discussion are, of course, welcome.