[WIP] Switch to unconstrained optimization followed by "rectifying" θ

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.35.2"
+version = "4.36.0"
 
 [deps]
 Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45"
@@ -14,11 +14,11 @@ JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 MixedModelsDatasets = "7e9fb7ac-9f67-43bf-b2c8-96ba0796cbb6"
-NLopt = "76087f3c-5699-56af-9a33-bf431cd00edd"
 PooledArrays = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+PRIMA = "0a7d04aa-8ac2-47b3-b7a7-9dbd6ad661ed"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
@@ -30,12 +30,6 @@ StructTypes = "856f2bd8-1eba-4b0a-8007-ebc267875bd4"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 TypedTables = "9d95f2ec-7b3d-5a63-8d20-e2491e220bb9"
 
-[weakdeps]
-PRIMA = "0a7d04aa-8ac2-47b3-b7a7-9dbd6ad661ed"
-
-[extensions]
-MixedModelsPRIMAExt = ["PRIMA"]
-
 [compat]
 Aqua = "0.8"
 Arrow = "1, 2"
@@ -51,7 +45,6 @@ JSON3 = "1"
 LinearAlgebra = "1"
 Markdown = "1"
 MixedModelsDatasets = "0.1"
-NLopt = "0.6, 1"
 PRIMA = "0.2"
 PooledArrays = "0.5, 1"
 PrecompileTools = "1"
@@ -77,10 +70,9 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
-PRIMA = "0a7d04aa-8ac2-47b3-b7a7-9dbd6ad661ed"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "DataFrames", "ExplicitImports", "InteractiveUtils", "PRIMA", "StableRNGs", "Suppressor", "Test"]
+test = ["Aqua", "DataFrames", "ExplicitImports", "InteractiveUtils", "StableRNGs", "Suppressor", "Test"]

issues/705/Project.toml

@@ -0,0 +1,2 @@
+[deps]
+MixedModels = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316"

src/MixedModels.jl

@@ -20,8 +20,8 @@ using LinearAlgebra: ldiv!, lmul!, logdet, mul!, norm, normalize, normalize!, qr
 using LinearAlgebra: rank, rdiv!, rmul!, svd, tril!
 using Markdown: Markdown
 using MixedModelsDatasets: dataset, datasets
-using NLopt: NLopt, Opt
 using PooledArrays: PooledArrays, PooledArray
+using PRIMA
 using PrecompileTools: PrecompileTools, @setup_workload, @compile_workload
 using ProgressMeter: ProgressMeter, Progress, ProgressUnknown, finish!, next!
 using Random: Random, AbstractRNG, randn!
@@ -128,10 +128,10 @@ export @formula,
     parametricbootstrap,
     pirls!,
     predict,
-    profile,
-    profileσ,
-    profilesigma,
-    profilevc,
+#    profile,
+#    profileσ,
+#    profilesigma,
+#    profilevc,
     pvalue,
     pwrss,
     ranef,
@@ -212,14 +212,13 @@ include("blockdescription.jl")
 include("grouping.jl")
 include("mimeshow.jl")
 include("serialization.jl")
-include("profile/profile.jl")
-include("nlopt.jl")
+#include("profile/profile.jl")
 include("prima.jl")
 
 
 # aliases with non-unicode function names
 const settheta! = setθ!
-const profilesigma = profileσ
+# const profilesigma = profileσ
 
 # COV_EXCL_START
 @setup_workload begin

src/linearmixedmodel.jl

@@ -496,10 +496,13 @@ function StatsAPI.fit!(
     end
 
     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
+    lb = optsum.lowerbd              # FIXME: this will not exist in the future.  Check for on diagonal instead
     for i in eachindex(xmin_)
         if iszero(lb[i]) && zero(T) < xmin_[i] < optsum.xtol_zero_abs
             xmin_[i] = zero(T)
@@ -517,7 +520,7 @@ function StatsAPI.fit!(
     # because we do that during the initial guess and rescale check
 
     ## ensure that the parameter values saved in m are xmin
-    objective!(m)(xmin)
+    objective!(m, xmin)
 
     optsum.final = xmin
     optsum.fmin = fmin
@@ -954,6 +957,36 @@ function PCA(m::LinearMixedModel; corr::Bool=true)
     return NamedTuple{_unique_fnames(m)}(PCA.(m.reterms, corr=corr))
 end
 
+"""
+    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)
+    foreach(m.λ) do λ
+        rectify!(λ)
+    end
+    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
+
 """
     reevaluateAend!(m::LinearMixedModel)
 

src/optsummary.jl

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