[WIP] Switch to unconstrained optimizers - step 1

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 <me@phillipalday.com>", "Douglas Bates <dmbates@gmail.com>"]
-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 @@
 
     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_)
@@ -735,7 +739,7 @@
     io::IO, ::MIME"text/plain", m::GeneralizedLinearMixedModel{T,D}
 ) where {T,D}
     if m.optsum.feval < 0
         @warn("Model has not been fit")
         return nothing
     end
     nAGQ = m.LMM.optsum.nAGQ
@@ -789,10 +793,10 @@
     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, :)

src/profile/thetapr.jl

@@ -28,12 +28,13 @@ function profileθj!(
 ) where {T}
     (; m, fwd, rev) = val
     optsum = m.optsum
-    (; final, fmin, lowerbd) = optsum
+    (; final, fmin) = optsum
     j = parsej(sym)
     θ = copy(final)
-    lbj = lowerbd[j]
     osj = optsum
     opt = Opt(osj)
+    pmj = m.parmap[j]
+    lbj = pmj[2] == pmj[3] ? zero(T) : T(-Inf)
     if length(θ) > 1      # set up the conditional optimization problem
         notj = deleteat!(collect(axes(final, 1)), j)
         osj = optsumj(optsum, j)

src/remat.jl

@@ -180,14 +180,12 @@ nθ(A::ReMat) = length(A.inds)
 """
     lowerbd{T}(A::ReMat{T})
 
-Return the vector of lower bounds on the parameters, `θ` associated with `A`
-
-These are the elements in the lower triangle of `A.λ` in column-major ordering.
-Diagonals have a lower bound of `0`.  Off-diagonals have a lower-bound of `-Inf`.
+Return the vector of lower bounds on the parameters, `θ` associated with `A`.  For unconstrained optimization these are all T(-Inf)
 """
 function lowerbd(A::ReMat{T}) where {T}
-    k = size(A.λ, 1)  # construct diagind(A.λ) by hand following #52115
-    return T[x ∈ range(1; step=k + 1, length=k) ? zero(T) : T(-Inf) for x in A.inds]
+    return fill!(similar(A.inds, T), -Inf)
+    #    k = size(A.λ, 1)  # construct diagind(A.λ) by hand following #52115
+    #    return T[x ∈ range(1; step=k + 1, length=k) ? T(-) : T(-Inf) for x in A.inds]
 end
 
 """