[WIP] Switch to unconstrained optimizers - step 1 (#840)

* Change lower bounds to -Inf, repair tests

* run JuliaFormatter on src dir

* Comparison values from AppleAccelerate

* [ci skip] allow PRIMA::newuoa optimizer

* [ci skip] don't pass lower bounds to optimizer

* Initial comparison of optimizers

* [ci skip] Adjust test targets and tolerances

* [ci skip] Adjust tolerance on test for x86_64, Linux, OpenBLAS

* Force ci

* Loosen tolerances (again)

* yet another tolerance adjustment in a test

* format

* Coverage and formatting

* format

* Add missing import

* Default optimizer to :LN_NEWUOA. Add timingtable function.

* Account for re-ordering of models.

* Default optimizer to :LN_NEWUOA, adjust tests, add timingtable

* Remove test of nonsensical model.

* Add atol to first(std(fmnc)) test.

* Add tolerance in prima test.

* More tolerance additions

* bump version to 5.0.0-DEV

* NEWS update

---------

Co-authored-by: Phillip Alday <[email protected]>

Author: dmbates

NEWS.md

@@ -1,3 +1,10 @@
+MixedModels v5.0.0 Release Notes
+==============================
+- Optimization is now performed _without constraints_. In a post-fitting step, the the Cholesky factor is canonicalized to have non-negative diagonal elements. [#840]
+- The default optimizer has changed to NLopt's implementation of NEWUOA where possible. NLopt's implementation fails on 1-dimensional problems, so in the case A single, scalar random effect, BOBYQA is used instead. In the future, the default optimizer backend will likely change to PRIMA and NLopt support will be moved to an extension. Blocking this change in backend is an issue with PRIMA.jl when running in VSCode's built-in REPL on Linux. [#840]
+- [BREAKING] Support for constrained optimization has been completely removed, i.e. the field `lowerbd` has been removed from `OptSummary`.
+
+
 MixedModels v4.38.0 Release Notes
 ==============================
 - Experimental support for evaluating `FiniteDiff.finite_difference_gradient` and `FiniteDiff.finite_difference_hessian of the objective of a fitted `LinearMixedModel`. [#842]
@@ -655,5 +662,6 @@ Package dependencies
 [#828]: https://github.com/JuliaStats/MixedModels.jl/issues/828
 [#829]: https://github.com/JuliaStats/MixedModels.jl/issues/829
 [#836]: https://github.com/JuliaStats/MixedModels.jl/issues/836
+[#840]: https://github.com/JuliaStats/MixedModels.jl/issues/840
 [#841]: https://github.com/JuliaStats/MixedModels.jl/issues/841
 [#842]: https://github.com/JuliaStats/MixedModels.jl/issues/842

Project.toml

@@ -1,7 +1,7 @@
 name = "MixedModels"
 uuid = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316"
 author = ["Phillip Alday <[email protected]>", "Douglas Bates <[email protected]>"]
-version = "4.38.0"
+version = "5.0.0-DEV"
 
 [deps]
 Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45"

ext/MixedModelsForwardDiffExt.jl

@@ -9,6 +9,7 @@ using MixedModels: AbstractReMat,
     kp1choose2,
     LD,
     lmulΛ!,
+    rankUpdate!,
     rmulΛ!,
     ssqdenom,
     UniformBlockDiagonal

ext/MixedModelsPRIMAExt.jl

@@ -25,6 +25,7 @@ end
 prima_optimizer!(::Val{:bobyqa}, args...; kwargs...) = PRIMA.bobyqa!(args...; kwargs...)
 prima_optimizer!(::Val{:cobyla}, args...; kwargs...) = PRIMA.cobyla!(args...; kwargs...)
 prima_optimizer!(::Val{:lincoa}, args...; kwargs...) = PRIMA.lincoa!(args...; kwargs...)
+prima_optimizer!(::Val{:newuoa}, args...; kwargs...) = PRIMA.newuoa!(args...; kwargs...)
 
 function MixedModels.optimize!(m::LinearMixedModel, ::PRIMABackend;
     progress::Bool=true, fitlog::Bool=false, kwargs...)
@@ -56,7 +57,8 @@ function MixedModels.optimize!(m::LinearMixedModel, ::PRIMABackend;
 
     maxfun = optsum.maxfeval > 0 ? optsum.maxfeval : 500 * length(optsum.initial)
     info = prima_optimizer!(Val(optsum.optimizer), obj, optsum.final;
-        xl=optsum.lowerbd, maxfun,
+        #        xl=optsum.lowerbd, 
+        maxfun,
         optsum.rhoend, optsum.rhobeg)
     ProgressMeter.finish!(prog)
     optsum.feval = info.nf

src/generalizedlinearmixedmodel.jl

@@ -296,6 +296,10 @@ function StatsAPI.fit!(
 
     xmin, fmin = optimize!(m; progress, fitlog, fast, verbose, nAGQ)
 
+    θopt = length(xmin) == length(θ) ? xmin : view(xmin, (length(β) + 1):lastindex(xmin))
+    rectify!(m.LMM)                  # flip signs of columns of m.λ elements with negative diagonal els
+    getθ!(θopt, m)                   # use the rectified values in xmin
+
     ## check if very small parameter values bounded below by zero can be set to zero
     xmin_ = copy(xmin)
     for i in eachindex(xmin_)
@@ -789,10 +793,10 @@ function unfit!(model::GeneralizedLinearMixedModel{T}) where {T}
     optsum = model.LMM.optsum
     # we need to reset optsum so that it
     # plays nice with the modifications fit!() does
-    optsum.lowerbd = mapfoldl(lowerbd, vcat, reterms)
+    optsum.lowerbd = mapfoldl(lowerbd, vcat, reterms)      # probably don't need this anymore - now trivial with all elements = -Inf
     # for variances (bounded at zero), we have ones, while
     # for everything else (bounded at -Inf), we have zeros
-    optsum.initial = map(T ∘ iszero, optsum.lowerbd)
+    optsum.initial = map(x -> T(x[2] == x[3]), model.LMM.parmap)
     optsum.final = copy(optsum.initial)
     optsum.xtol_abs = fill!(copy(optsum.initial), 1.0e-10)
     optsum.initial_step = T[]

src/linearmixedmodel.jl

@@ -497,11 +497,14 @@ function StatsAPI.fit!(
 
     xmin, fmin = optimize!(m; progress, fitlog)
 
+    setθ!(m, xmin)                   # ensure that the parameters saved in m are xmin
+    rectify!(m)                      # flip signs of columns of m.λ elements with negative diagonal els
+    getθ!(xmin, m)                   # use the rectified values as xmin
+
     ## check if small non-negative parameter values can be set to zero
     xmin_ = copy(xmin)
-    lb = optsum.lowerbd
-    for i in eachindex(xmin_)
-        if iszero(lb[i]) && zero(T) < xmin_[i] < optsum.xtol_zero_abs
+    for (i, pm) in enumerate(m.parmap)
+        if pm[2] == pm[3] && zero(T) < xmin_[i] < optsum.xtol_zero_abs
             xmin_[i] = zero(T)
         end
     end
@@ -681,6 +684,19 @@ function Base.getproperty(m::LinearMixedModel{T}, s::Symbol) where {T}
     end
 end
 
+"""
+    _set_init(m::LinearMixedModel)
+
+Set each element of m.optsum.initial to 1.0 for diagonal and 0.0 for off-diagonal
+"""
+function _set_init!(m::LinearMixedModel)
+    init = m.optsum.initial
+    for (i, pm) in enumerate(m.parmap)
+        init[i] = pm[2] == pm[3]
+    end
+    return m
+end
+
 StatsAPI.islinear(m::LinearMixedModel) = true
 
 """
@@ -927,6 +943,34 @@ end
 
 LinearAlgebra.rank(m::LinearMixedModel) = m.feterm.rank
 
+"""
+    rectify!(m::LinearMixedModel)
+
+For each element of m.λ check for negative values on the diagonal and flip the signs of the entire column when any are present.
+
+This provides a canonical converged value of θ.  We use unconstrained optimization followed by this reassignment to avoid the
+hassle of constrained optimization.
+"""
+function rectify!(m::LinearMixedModel)
+    rectify!.(m.λ)
+    return m
+end
+
+function rectify!(λ::LowerTriangular)
+    for (j, c) in enumerate(eachcol(λ.data))
+        if c[j] < 0
+            c .*= -1
+        end
+    end
+    return λ
+end
+
+function rectify!(λ::Diagonal)
+    d = λ.diag
+    map!(abs, d, d)
+    return λ
+end
+
 """
     rePCA(m::LinearMixedModel; corr::Bool=true)
 
@@ -939,7 +983,7 @@ principal component, the first two principal components, etc.  The last element
 always 1.0 representing the complete proportion of the variance.
 """
 function rePCA(m::LinearMixedModel; corr::Bool=true)
-    pca = PCA.(m.reterms, corr=corr)
+    pca = PCA.(m.reterms; corr=corr)
     return NamedTuple{_unique_fnames(m)}(getproperty.(pca, :cumvar))
 end
 
@@ -951,7 +995,7 @@ covariance matrices or correlation matrices when `corr` is `true`.
 """
 
 function PCA(m::LinearMixedModel; corr::Bool=true)
-    return NamedTuple{_unique_fnames(m)}(PCA.(m.reterms, corr=corr))
+    return NamedTuple{_unique_fnames(m)}(PCA.(m.reterms; corr=corr))
 end
 
 """
@@ -1256,9 +1300,8 @@ function unfit!(model::LinearMixedModel{T}) where {T}
     optsum = model.optsum
     optsum.feval = -1
     optsum.initial_step = T[]
-    # for variances (bounded at zero), we have ones, while
-    # for everything else (bounded at -Inf), we have zeros
-    map!(T ∘ iszero, optsum.initial, optsum.lowerbd)
+    # initialize elements on the diagonal of Λ to one(T), off-diagonals to zero(T)
+    _set_init!(model)
     copyto!(optsum.final, optsum.initial)
     reevaluateAend!(model)
 

src/mixedmodel.jl

@@ -68,8 +68,11 @@ Equality comparisons are used b/c small non-negative θ values are replaced by 0
     For `GeneralizedLinearMixedModel`, the entire parameter vector (including
     β in the case `fast=false`) must be specified if the default is not used.
 """
-function issingular(m::MixedModel, θ=m.θ; atol::Real=0, rtol::Real=atol > 0 ? 0 : √eps())
-    return _issingular(m.lowerbd, θ; atol, rtol)
+function issingular(
+    m::MixedModel{T}, θ=m.θ; atol::Real=0, rtol::Real=atol > 0 ? 0 : √eps()
+) where {T}
+    lb = [(pm[2] == pm[3]) ? zero(T) : T(-Inf) for pm in m.parmap]
+    return _issingular(lb, θ; atol, rtol)
 end
 
 function _issingular(v, w; atol, rtol)

src/nlopt.jl

@@ -82,15 +82,19 @@ end
 
 function NLopt.Opt(optsum::OptSummary)
     lb = optsum.lowerbd
+    n = length(lb)
 
-    opt = NLopt.Opt(optsum.optimizer, length(lb))
+    if optsum.optimizer == :LN_NEWUOA && isone(n) # :LN_NEWUOA doesn't allow n == 1
+        optsum.optimizer = :LN_BOBYQA
+    end
+    opt = NLopt.Opt(optsum.optimizer, n)
     NLopt.ftol_rel!(opt, optsum.ftol_rel) # relative criterion on objective
     NLopt.ftol_abs!(opt, optsum.ftol_abs) # absolute criterion on objective
     NLopt.xtol_rel!(opt, optsum.xtol_rel) # relative criterion on parameter values
-    if length(optsum.xtol_abs) == length(lb)  # not true for fast=false optimization in GLMM
+    if length(optsum.xtol_abs) == n  # not true for fast=false optimization in GLMM
         NLopt.xtol_abs!(opt, optsum.xtol_abs) # absolute criterion on parameter values
     end
-    NLopt.lower_bounds!(opt, lb)
+    #    NLopt.lower_bounds!(opt, lb)   # use unconstrained optimization even for :LN_BOBYQA
     NLopt.maxeval!(opt, optsum.maxfeval)
     NLopt.maxtime!(opt, optsum.maxtime)
     if isempty(optsum.initial_step)

src/optsummary.jl

@@ -60,7 +60,7 @@ Base.@kwdef mutable struct OptSummary{T<:AbstractFloat}
     ftol_zero_abs::T = eltype(initial)(1.e-5)
     maxfeval::Int = -1
 
-    optimizer::Symbol = :LN_BOBYQA
+    optimizer::Symbol = :LN_NEWUOA    # switched to :LN_BOBYQA for one-dimensional optimizations
     backend::Symbol = :nlopt
 
     # the @kwdef macro isn't quite smart enough for us to use the type parameter
@@ -85,7 +85,7 @@ end
 function OptSummary(
     initial::Vector{T},
     lowerbd::Vector{S},
-    optimizer::Symbol=:LN_BOBYQA; kwargs...,
+    optimizer::Symbol=:LN_NEWUOA; kwargs...,
 ) where {T<:AbstractFloat,S<:AbstractFloat}
     TS = promote_type(T, S)
     return OptSummary{TS}(; initial, lowerbd, optimizer, kwargs...)

src/pca.jl

@@ -81,7 +81,7 @@ function Base.show(
         # only display the lower triangle of symmetric matrix
         if pca.rnames !== missing
             n = length(pca.rnames)
-            cv = string.(round.(pca.covcor, digits=ndigitsmat))
+            cv = string.(round.(pca.covcor; digits=ndigitsmat))
             dotpad = lpad(".", div(maximum(length, cv), 2))
             for i in 1:n, j in (i + 1):n
                 cv[i, j] = dotpad
@@ -97,7 +97,7 @@ function Base.show(
             # if there are no names, then we cheat and use the print method
             # for LowerTriangular, which automatically covers the . in the
             # upper triangle
-            printmat = round.(LowerTriangular(pca.covcor), digits=ndigitsmat)
+            printmat = round.(LowerTriangular(pca.covcor); digits=ndigitsmat)
         end
 
         Base.print_matrix(io, printmat)
@@ -106,21 +106,21 @@ function Base.show(
     if stddevs
         println(io, "\nStandard deviations:")
         sv = pca.sv
-        show(io, round.(sv.S, digits=ndigitsvec))
+        show(io, round.(sv.S; digits=ndigitsvec))
         println(io)
     end
     if variances
         println(io, "\nVariances:")
         vv = abs2.(sv.S)
-        show(io, round.(vv, digits=ndigitsvec))
+        show(io, round.(vv; digits=ndigitsvec))
         println(io)
     end
     println(io, "\nNormalized cumulative variances:")
-    show(io, round.(pca.cumvar, digits=ndigitscum))
+    show(io, round.(pca.cumvar; digits=ndigitscum))
     println(io)
     if loadings
         println(io, "\nComponent loadings")
-        printmat = round.(pca.loadings, digits=ndigitsmat)
+        printmat = round.(pca.loadings; digits=ndigitsmat)
         if pca.rnames !== missing
             pclabs = [Text(""); Text.("PC$i" for i in 1:length(pca.rnames))]
             pclabs = reshape(pclabs, 1, :)