JuliaStats
Different results with Julia 1.11 on macOS
The following example gives a different result with Julia 1.11.5 compared with Julia 1.10.5
using DataFrames, CSV, MixedModels
df = CSV.read("data.csv", DataFrame)
_lmm = LinearMixedModel(
@formula(endpoint ~ 1 + formulation + sequence + period + (1 | id)),
df;
contrasts = Dict(:period => DummyCoding()),
)
_lmm.optsum.optimizer = :LN_COBYLA
fit!(_lmm; REML = true)
This doesn't seem to happen on Linux so it might be related to minor rounding differences. If I use BOBYQA then I'm getting the better solution. We are using COBYLA because of https://github.com/JuliaStats/MixedModels.jl/issues/705.
Julia 1.10.5
Linear mixed model fit by REML
endpoint ~ 1 + formulation + sequence + period + (1 | id)
REML criterion at convergence: 1143.3630490413652
Variance components:
Column Variance Std.Dev.
id (Intercept) 0.842408 0.917828
Residual 2.118190 1.455400
Number of obs: 298; levels of grouping factors: 77
Fixed-effects parameters:
───────────────────────────────────────────────────────
Coef. Std. Error z Pr(>|z|)
───────────────────────────────────────────────────────
(Intercept) -0.313961 0.254152 -1.24 0.2167
formulation: T 0.177572 0.168966 1.05 0.2933
sequence: TRTR 0.330212 0.269183 1.23 0.2199
period: 2 -1.39972 0.235472 -5.94 <1e-08
period: 3 0.947798 0.241498 3.92 <1e-04
period: 4 -0.0445583 0.236472 -0.19 0.8505
───────────────────────────────────────────────────────
julia> _lmm.optsum
...
Function evaluations: 28
xtol_zero_abs: 0.001
ftol_zero_abs: 1.0e-5
Final parameter vector: [0.6306362152099609]
Final objective value: 1143.3630490413652
Return code: FTOL_REACHED
julia> _lmm.optsum.fitlog
28-element Vector{Tuple{Vector{Float64}, Float64}}:
([1.0], 1153.9741825642604)
([1.75], 1197.7827103789102)
([0.2500000000000001], 1159.431337870577)
([0.625], 1143.3663690793285)
([0.24999999999999994], 1159.431337870577)
([0.8125], 1146.367205373148)
([0.53125], 1144.4529917605419)
([0.671875], 1143.5350688209273)
([0.6015625], 1143.4526930723844)
([0.63671875], 1143.3668820226)
([0.619140625], 1143.376909665865)
⋮
([0.6304931640625], 1143.3630512130903)
([0.630584716796875], 1143.3630493324415)
([0.6307220458984375], 1143.3630497829142)
([0.6306533813476562], 1143.3630490670132)
([0.6306304931640625], 1143.363049046446)
([0.6306076049804688], 1143.3630491349204)
([0.6306190490722656], 1143.3630490770527)
([0.6306362152099609], 1143.3630490413652)
([0.6306419372558594], 1143.3630490430996)
([0.6306390762329102], 1143.3630490413805)
Julia 1.11.5
Linear mixed model fit by REML
endpoint ~ 1 + formulation + sequence + period + (1 | id)
REML criterion at convergence: 1143.3663690793285
Variance components:
Column Variance Std.Dev.
id (Intercept) 0.829727 0.910894
Residual 2.124101 1.457430
Number of obs: 298; levels of grouping factors: 77
Fixed-effects parameters:
───────────────────────────────────────────────────────
Coef. Std. Error z Pr(>|z|)
───────────────────────────────────────────────────────
(Intercept) -0.313958 0.253726 -1.24 0.2159
formulation: T 0.177514 0.169199 1.05 0.2941
sequence: TRTR 0.330265 0.268101 1.23 0.2180
period: 2 -1.39962 0.235799 -5.94 <1e-08
period: 3 0.947926 0.241825 3.92 <1e-04
period: 4 -0.0446212 0.236798 -0.19 0.8505
───────────────────────────────────────────────────────
julia> _lmm.optsum
...
Function evaluations: 5
xtol_zero_abs: 0.001
ftol_zero_abs: 1.0e-5
Final parameter vector: [0.625]
Final objective value: 1143.3663690793285
Return code: FTOL_REACHED
julia> _lmm.optsum.fitlog
5-element Vector{Tuple{Vector{Float64}, Float64}}:
([1.0], 1153.9741825642607)
([1.75], 1197.7827103789102)
([0.25], 1159.4313378705772)
([0.625], 1143.3663690793285)
([0.25000000000000006], 1159.431337870577)
Versioninfo
julia> versioninfo()
Julia Version 1.11.5
Commit 760b2e5b739 (2025-04-14 06:53 UTC)
Build Info:
Official https://julialang.org/ release
Platform Info:
OS: macOS (arm64-apple-darwin24.0.0)
CPU: 8 × Apple M1 Pro
WORD_SIZE: 64
LLVM: libLLVM-16.0.6 (ORCJIT, apple-m1)
Threads: 1 default, 0 interactive, 1 GC (on 6 virtual cores)
Update: I've added the fitlog values. I think it looks like an issue in nlops's COBYLA. Why would it stop with FTOL_REACHED after five iterations. Non of the values are that close.
@andreasnoack Have you tried the PRIMA versions of BOBYQA and COBYLA? There is an extension to MixedModels.jl that is available if you add PRIMA first in your session. It provides a prfit! function. I'll see if I can reproduce these results with PRIMA. The author of libprima said that he made some code corrections when he transcribed Powell's Fortran77 code to a more modern version of Fortran.
Indeed, I started playing around with PRIMA. It looks like the PRIMA COBYLA works here while PRIMA's BOBYQA seems to have inherited the issue from Powell/nlopt, see https://github.com/JuliaStats/MixedModels.jl/issues/705#issuecomment-2980446198. I should also mention that I hit some warnings that the maximum number of function evaluations had been reached with PRIMA's COBYLA but I simply increased optsum.maxfeval and that was it.
I agree we need to look into this in greater detail. The issue @ajinkya-k mentioned is that editors trying to reproduce results from a paper we submitted got substantially different results. However, we don't yet know what versions of Julia, MixedModels, etc. they were using.
I did try switching from LN_BOBYQA to LN_COBYLA on that model fit in Julia 1.10.9 and got essentially the same result but at a cost of 2446 function evaluations as opposed to 170 for BOBYQA. See enclosed
Another issue for us is the use of accelerated BLAS, such as AppleAccelerate or MKL. When you have two or more grouping factors for the random effects and large numbers of levels for them the second diagonal block in the blocked Cholesky factor usually ends up being large and dense so accelerated BLAS have a big impact. Because they often re-order the arithmetic operations the value of the objective can depend on details of the processor, resulting in slightly different optima.
Our example is described in https://arxiv.org/abs/2505.11674
I tried the PRIMA versions of BOBYQA and COBYLA under Julia 1.10.9 with essentially the same results - convergence to an objective of 237648.6016 after 164 evaluations for BOBYQA and 1402 evaluations for COBYLA.
I tried to recreate the premature termination with Julia-1.11.5 on an M-series MacOS laptop and was unable to
julia> using DataFrames, CSV, MixedModels, CairoMakie
julia> df = CSV.read("issues/833/data.csv", DataFrame)
298×6 DataFrame
Row │ id id_within_sequence sequence period formulation endpoint
│ Int64 Int64 String7 Int64 String1 Float64
─────┼──────────────────────────────────────────────────────────────────────
1 │ 1 1 RTRT 1 R 0.237283
2 │ 1 1 RTRT 2 T -1.09763
3 │ 1 1 RTRT 3 R 2.81768
4 │ 1 1 RTRT 4 T 1.0022
⋮ │ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮
296 │ 78 38 RTRT 2 T -1.82628
297 │ 78 38 RTRT 3 R -2.00504
298 │ 78 38 RTRT 4 T 1.10657
291 rows omitted
julia> _lmm = LinearMixedModel(
@formula(endpoint ~ 1 + formulation + sequence + period + (1 | id)),
df;
contrasts = Dict(:period => DummyCoding()),
)
┌ Warning: Model has not been fit
└ @ MixedModels ~/.julia/packages/MixedModels/TQ7b4/src/linearmixedmodel.jl:1056
julia> _lmm.optsum.optimizer = :LN_COBYLA
:LN_COBYLA
julia> fit!(_lmm; REML = true, fitlog=true)
Linear mixed model fit by REML
endpoint ~ 1 + formulation + sequence + period + (1 | id)
REML criterion at convergence: 1143.36304904116
Variance components:
Column VarianceStd.Dev.
id (Intercept) 0.84241 0.91783
Residual 2.11819 1.45540
Number of obs: 298; levels of grouping factors: 77
Fixed-effects parameters:
───────────────────────────────────────────────────────
Coef. Std. Error z Pr(>|z|)
───────────────────────────────────────────────────────
(Intercept) -0.313961 0.254152 -1.24 0.2167
formulation: T 0.177572 0.168966 1.05 0.2933
sequence: TRTR 0.330212 0.269183 1.23 0.2199
period: 2 -1.39972 0.235472 -5.94 <1e-08
period: 3 0.947798 0.241498 3.92 <1e-04
period: 4 -0.0445583 0.236472 -0.19 0.8505
───────────────────────────────────────────────────────
julia> lines(0.4:0.01:1.0, x -> objective(updateL!(setθ!(_lmm, (x,)))))
julia> versioninfo()
Julia Version 1.11.5
Commit 760b2e5b739 (2025-04-14 06:53 UTC)
Build Info:
Official https://julialang.org/ release
Platform Info:
OS: macOS (arm64-apple-darwin24.0.0)
CPU: 8 × Apple M1 Pro
WORD_SIZE: 64
LLVM: libLLVM-16.0.6 (ORCJIT, apple-m1)
Threads: 6 default, 0 interactive, 3 GC (on 6 virtual cores)
julia> _lmm.optsum.fitlog
32-element Vector{Tuple{Vector{Float64}, Float64}}:
([1.0], 1153.9741825642604)
([1.75], 1197.7827103789102)
([0.2500000000000001], 1159.431337870577)
([0.625], 1143.3663690793285)
([0.24999999999999994], 1159.431337870577)
([0.8125], 1146.3672053731477)
⋮
([0.6306362152099609], 1143.3630490413652)
([0.6306347846984863], 1143.3630490419964)
([0.6306369304656982], 1143.3630490412095)
([0.6306376457214355], 1143.36304904116)
([0.6306383609771729], 1143.3630490412172)
julia> _lmm.optsum
Initial parameter vector: [1.0]
Initial objective value: 1153.9741825642604
Backend: nlopt
Optimizer: LN_COBYLA
Lower bounds: [0.0]
ftol_rel: 1.0e-12
ftol_abs: 1.0e-8
xtol_rel: 0.0
xtol_abs: [1.0e-10]
initial_step: [0.75]
maxfeval: -1
maxtime: -1.0
Function evaluations: 32
xtol_zero_abs: 0.001
ftol_zero_abs: 1.0e-5
Final parameter vector: [0.6306376457214355]
Final objective value: 1143.36304904116
Return code: FTOL_REACHED
I added a plot of the objective, which is very close to quadratic, so this should not be a tricky optimization problem. Paradoxically one-dimensional optimization can be difficult for general optimization algorithms. BOBYQA is generally good on one-dimensional problems, which is one of the reasons that we use it as a default.
I have created a db/issue833 branch with the data and the script I used in an issues/833 directory.
With MixedModels v5.0.0 and later the default optimizer, which is currently :LN_BOBYQA (without lower bounds) for 1-dimentional problems, seems to work fine on Mac and Linux.
@andreasnoack Do you think this issue can be closed now?
Indeed. I've tested with the new version and this issue seems to be resolved.