StructuralEquationModels · Maximilian-Stefan-Ernst · Mar 22, 2025 · Mar 22, 2025 · Mar 22, 2025 · Mar 22, 2025
diff --git a/Project.toml b/Project.toml
@@ -8,6 +8,7 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
+InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LazyArtifacts = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
 LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -30,6 +31,7 @@ StenoGraphs = "0.2 - 0.3, 0.4.1 - 0.5"
 DataFrames = "1"
 Distributions = "0.25"
 FiniteDiff = "2"
+InteractiveUtils = "1.11.0"
 LineSearches = "7"
 NLSolversBase = "7"
 NLopt = "0.6, 1"

diff --git a/docs/make.jl b/docs/make.jl
@@ -1,6 +1,12 @@
 using Documenter, StructuralEquationModels
+using NLopt, ProximalAlgorithms
 
 makedocs(
+    modules=[
+        StructuralEquationModels,
+        Base.get_extension(StructuralEquationModels, :SEMNLOptExt),
+        Base.get_extension(StructuralEquationModels, :SEMProximalOptExt)
+        ],
     sitename = "StructuralEquationModels.jl",
     pages = [
         "index.md",
@@ -60,6 +66,7 @@ makedocs(
         collapselevel = 1,
     ),
     doctest = false,
+    checkdocs = :none,
 )
 
 # doctest(StructuralEquationModels, fix=true)

diff --git a/docs/src/developer/optimizer.md b/docs/src/developer/optimizer.md
@@ -30,7 +30,6 @@ update_observed(optimizer::SemOptimizerName, observed::SemObserved; kwargs...) =
 ### additional methods
 ############################################################################################
 
-algorithm(optimizer::SemOptimizerName) = optimizer.algorithm
 options(optimizer::SemOptimizerName) = optimizer.options
 ```
 
@@ -68,7 +67,7 @@ The method has to return a `SemFit` object that consists of the minimum of the o
 In addition, you might want to provide methods to access properties of your optimization result:
 
 ```julia
-optimizer(res::MyOptimizationResult) = ...
+algorithm_name(res::MyOptimizationResult) = ...
 n_iterations(res::MyOptimizationResult) = ...
 convergence(res::MyOptimizationResult) = ...
 ```
diff --git a/docs/src/performance/simulation.md b/docs/src/performance/simulation.md
@@ -67,7 +67,7 @@ For example,
 
 new_observed = SemObservedData(;data = data_2, specification = partable)
 
-my_optimizer = SemOptimizerOptim()
+my_optimizer = SemOptimizer()
 
 new_optimizer = update_observed(my_optimizer, new_observed)
 ```

diff --git a/docs/src/tutorials/backends/nlopt.md b/docs/src/tutorials/backends/nlopt.md
@@ -1,47 +1,36 @@
 # Using NLopt.jl
 
-[`SemOptimizerNLopt`](@ref) implements the connection to `NLopt.jl`.
-It is only available if the `NLopt` package is loaded alongside `StructuralEquationModels.jl` in the running Julia session.
-It takes a bunch of arguments:
+When [`NLopt.jl`](https://github.com/jump-dev/NLopt.jl) is loaded in the running Julia session,
+it could be used by the [`SemOptimizer`](@ref) by specifying `engine = :NLopt`
+(see ...).
+Among other things, `NLopt` enables constrained optimization of the SEM models, which is
+explained in the [Constrained optimization](@ref) section.
 
-```julia
-    •  algorithm: optimization algorithm
-
-    •  options::Dict{Symbol, Any}: options for the optimization algorithm
-
-    •  local_algorithm: local optimization algorithm
-
-    •  local_options::Dict{Symbol, Any}: options for the local optimization algorithm
-
-    •  equality_constraints::Vector{NLoptConstraint}: vector of equality constraints
-
-    •  inequality_constraints::Vector{NLoptConstraint}: vector of inequality constraints
-```
-Constraints are explained in the section on [Constrained optimization](@ref).
-
-The defaults are LBFGS as the optimization algorithm and the standard options from `NLopt.jl`.
-We can choose something different:
+We can override the default *NLopt* algorithm (LFBGS) and instead use
+the *augmented lagrangian* method with LBFGS as the *local* optimization algorithm,
+stop at a maximum of 200 evaluations and use a relative tolerance of
+the objective value of `1e-6` as the stopping criterion for the local algorithm:
 
 ```julia
 using NLopt
 
-my_optimizer = SemOptimizerNLopt(;
+my_optimizer = SemOptimizer(;
+    engine = :NLopt,
     algorithm = :AUGLAG,
     options = Dict(:maxeval => 200),
     local_algorithm = :LD_LBFGS,
     local_options = Dict(:ftol_rel => 1e-6)
 )
 ```
 
-This uses an augmented lagrangian method with LBFGS as the local optimization algorithm, stops at a maximum of 200 evaluations and uses a relative tolerance of the objective value of `1e-6` as the stopping criterion for the local algorithm.
-
 To see how to use the optimizer to actually fit a model now, check out the [Model fitting](@ref) section.
 
-In the NLopt docs, you can find explanations about the different [algorithms](https://nlopt.readthedocs.io/en/latest/NLopt_Algorithms/) and a [tutorial](https://nlopt.readthedocs.io/en/latest/NLopt_Introduction/) that also explains the different options.
+In the *NLopt* docs, you can find details about the [optimization algorithms](https://nlopt.readthedocs.io/en/latest/NLopt_Algorithms/),
+and the [tutorial](https://nlopt.readthedocs.io/en/latest/NLopt_Introduction/) that demonstrates how to tweak their behavior.
 
 To choose an algorithm, just pass its name without the 'NLOPT\_' prefix (for example, 'NLOPT\_LD\_SLSQP' can be used by passing `algorithm = :LD_SLSQP`).
 
-The README of the [julia package](https://github.com/JuliaOpt/NLopt.jl) may also be helpful, and provides a list of options:
+The README of the [*NLopt.jl*](https://github.com/JuliaOpt/NLopt.jl) may also be helpful, and provides a list of options:
 
  - `algorithm`
  - `stopval`

diff --git a/docs/src/tutorials/backends/optim.md b/docs/src/tutorials/backends/optim.md
@@ -1,23 +1,23 @@
 # Using Optim.jl
 
-[`SemOptimizerOptim`](@ref) implements the connection to `Optim.jl`.
-It takes two arguments, `algorithm` and `options`.
-The defaults are LBFGS as the optimization algorithm and the standard options from `Optim.jl`.
-We can load the `Optim` and `LineSearches` packages to choose something different:
+[Optim.jl](https://github.com/JuliaNLSolvers/Optim.jl) is the default optimization engine of *SEM.jl*,
+see ... for a full list of its parameters.
+It defaults to the LBFGS optimization, but we can load the `Optim` and `LineSearches` packages
+and specify BFGS (!not L-BFGS) with a back-tracking linesearch and Hager-Zhang initial step length guess:
 
 ```julia
 using Optim, LineSearches
 
-my_optimizer = SemOptimizerOptim(
+my_optimizer = SemOptimizer(
     algorithm = BFGS(
-        linesearch = BackTracking(order=3), 
+        linesearch = BackTracking(order=3),
         alphaguess = InitialHagerZhang()
-        ),
-    options = Optim.Options(show_trace = true) 
-    )
+    ),
+    options = Optim.Options(show_trace = true)
+)
 ```
 
-This optimizer will use BFGS (!not L-BFGS) with a back tracking linesearch and a certain initial step length guess. Also, the trace of the optimization will be printed to the console.
+Note that we used `options` to print the optimization progress to the console.
 
 To see how to use the optimizer to actually fit a model now, check out the [Model fitting](@ref) section.
 

diff --git a/docs/src/tutorials/concept.md b/docs/src/tutorials/concept.md
@@ -21,13 +21,13 @@ So everything that can be used as the 'observed' part has to be of type `SemObse
 
 Here is an overview on the available building blocks:
 
-|[`SemObserved`](@ref)            | [`SemImplied`](@ref)  | [`SemLossFunction`](@ref) | [`SemOptimizer`](@ref)        |
-|---------------------------------|-----------------------|---------------------------|-------------------------------|
-| [`SemObservedData`](@ref)       | [`RAM`](@ref)         | [`SemML`](@ref)           | [`SemOptimizerOptim`](@ref)   |
-| [`SemObservedCovariance`](@ref) | [`RAMSymbolic`](@ref) | [`SemWLS`](@ref)          | [`SemOptimizerNLopt`](@ref)   |
-| [`SemObservedMissing`](@ref)    | [`ImpliedEmpty`](@ref)| [`SemFIML`](@ref)         |                               |
-|                                 |                       | [`SemRidge`](@ref)        |                               |
-|                                 |                       | [`SemConstant`](@ref)     |                               |
+|[`SemObserved`](@ref)            | [`SemImplied`](@ref)  | [`SemLossFunction`](@ref) | [`SemOptimizer`](@ref)     |
+|---------------------------------|-----------------------|---------------------------|----------------------------|
+| [`SemObservedData`](@ref)       | [`RAM`](@ref)         | [`SemML`](@ref)           | :Optim                     |
+| [`SemObservedCovariance`](@ref) | [`RAMSymbolic`](@ref) | [`SemWLS`](@ref)          | :NLopt                     |
+| [`SemObservedMissing`](@ref)    | [`ImpliedEmpty`](@ref)| [`SemFIML`](@ref)         | :Proximal                  |
+|                                 |                       | [`SemRidge`](@ref)        |                            |
+|                                 |                       | [`SemConstant`](@ref)     |                            |
 
 The rest of this page explains the building blocks for each part. First, we explain every part and give an overview on the different options that are available. After that, the [API - model parts](@ref) section serves as a reference for detailed explanations about the different options.
 (How to stick them together to a final model is explained in the section on [Model Construction](@ref).)
@@ -52,7 +52,7 @@ Available loss functions are
 ## The optimizer part aka `SemOptimizer`
 The optimizer part of a model connects to the numerical optimization backend used to fit the model.
 It can be used to control options like the optimization algorithm, linesearch, stopping criteria, etc.
-There are currently three available backends, [`SemOptimizerOptim`](@ref) connecting to the [Optim.jl](https://github.com/JuliaNLSolvers/Optim.jl) backend, [`SemOptimizerNLopt`](@ref) connecting to the [NLopt.jl](https://github.com/JuliaOpt/NLopt.jl) backend and [`SemOptimizerProximal`](@ref) connecting to [ProximalAlgorithms.jl](https://github.com/JuliaFirstOrder/ProximalAlgorithms.jl).
+There are currently three available engines (i.e., backends used to carry out the numerical optimization), `:Optim` connecting to the [Optim.jl](https://github.com/JuliaNLSolvers/Optim.jl) backend, `:NLopt` connecting to the [NLopt.jl](https://github.com/JuliaOpt/NLopt.jl) backend and `:Proximal` connecting to [ProximalAlgorithms.jl](https://github.com/JuliaFirstOrder/ProximalAlgorithms.jl).
 For more information about the available options see also the tutorials about [Using Optim.jl](@ref) and [Using NLopt.jl](@ref), as well as [Constrained optimization](@ref) and [Regularization](@ref) .
 
 # What to do next
@@ -102,7 +102,13 @@ SemConstant
 
 ```@docs
 SemOptimizer
-SemOptimizerOptim
-SemOptimizerNLopt
-SemOptimizerProximal
+```
+
+A reference: [NLopt engine](@ref SEMNLOptExt.SemOptimizerNLopt)
+
+```@autodocs
+Modules = [
+    Base.get_extension(StructuralEquationModels, :SEMNLOptExt),
+    Base.get_extension(StructuralEquationModels, :SEMProximalOptExt)]
+Order = [:type, :function]
 ```
diff --git a/docs/src/tutorials/constraints/constraints.md b/docs/src/tutorials/constraints/constraints.md
@@ -1,10 +1,15 @@
 # Constrained optimization
 
-## Using the NLopt backend
+*SEM.jl* allows to fit models with additional constraints imposed on the parameters.
+
+## Using the NLopt engine
+
+*NLopt.jl* is one of *SEM.jl* optimization engines that supports constrained optimization.
+In the example below we show how to specify constraints for the *SEM* model when using *NLopt*.
 
 ### Define an example model
 
-Let's revisit our model from [A first model](@ref):
+Let's revisit our model from [A first model](@ref) and fit it first without constraints:
 
 ```@example constraints
 using StructuralEquationModels
@@ -57,39 +62,40 @@ details(partable)
 
 ### Define the constraints
 
-Let's introduce some constraints:
+Let's introduce some constraints (they are not based on any real properties of the underlying study and serve only as an example):
 1. **Equality constraint**: The covariances `y3 ↔ y7` and `y8 ↔ y4` should sum up to `1`.
 2. **Inequality constraint**: The difference between the loadings `dem60 → y2` and `dem60 → y3` should be smaller than `0.1`
 3. **Bound constraint**: The directed effect from  `ind60 → dem65` should be smaller than `0.5`
 
-(Of course those constaints only serve an illustratory purpose.)
-
-We first need to get the indices of the respective parameters that are invoved in the constraints.
-We can look up their labels in the output above, and retrieve their indices as
+Since *NLopt* does not have access to the SEM parameter names, its constaints are defined on the vector of all SEM parameters.
+We have to look up the indices of the parameters involved in the constraints to construct the respective functions.
 
 ```@example constraints
 parind = param_indices(model)
 parind[:y3y7] # 29
 ```
 
-The bound constraint is easy to specify: Just give a vector of upper or lower bounds that contains the bound for each parameter. In our example, only the parameter labeled `:λₗ` has an upper bound, and the number of total parameters is `n_par(model) = 31`, so we define
+The bound constraint is easy to specify: just give a vector of upper or lower bounds for each parameter.
+In our example, only the parameter labeled `:λₗ` has an upper bound, and the number of total parameters is `n_par(model) = 31`, so
 
 ```@example constraints
 upper_bounds = fill(Inf, 31)
 upper_bounds[parind[:λₗ]] = 0.5
 ```
 
-The equailty and inequality constraints have to be reformulated to be of the form `x = 0` or `x ≤ 0`:
-1. `y3 ↔ y7 + y8 ↔ y4 - 1 = 0`
-2. `dem60 → y2 - dem60 → y3 - 0.1 ≤ 0`
+The equailty and inequality constraints have to be reformulated in the `f(θ) = 0` or `f(θ) ≤ 0` form,
+where `θ` is the vector of SEM parameters:
+1. `f(θ) = 0`, where `f(θ) = y3 ↔ y7 + y8 ↔ y4 - 1`
+2. `g(θ) ≤ 0`, where `g(θ) = dem60 → y2 - dem60 → y3 - 0.1`
 
-Now they can be defined as functions of the parameter vector:
+If the optimization algorithm needs gradients, it will pass the `gradient` vector that is of the same size as the parameters,
+and the constraint function has to calculate the gradient in-place.
 
 ```@example constraints
 parind[:y3y7] # 29
 parind[:y8y4] # 30
 # θ[29] + θ[30] - 1 = 0.0
-function eq_constraint(θ, gradient)
+function f(θ, gradient)
     if length(gradient) > 0
         gradient .= 0.0
         gradient[29] = 1.0
@@ -101,7 +107,7 @@ end
 parind[:λ₂] # 3
 parind[:λ₃] # 4
 # θ[3] - θ[4] - 0.1 ≤ 0
-function ineq_constraint(θ, gradient)
+function g(θ, gradient)
     if length(gradient) > 0
         gradient .= 0.0
         gradient[3] = 1.0
@@ -111,49 +117,43 @@ function ineq_constraint(θ, gradient)
 end
 ```
 
-If the algorithm needs gradients at an iteration, it will pass the vector `gradient` that is of the same size as the parameters.
-With `if length(gradient) > 0` we check if the algorithm needs gradients, and if it does, we fill the `gradient` vector with the gradients
-of the constraint w.r.t. the parameters.
-
-In NLopt, vector-valued constraints are also possible, but we refer to the documentation for that.
+In *NLopt*, vector-valued constraints are also possible, but we refer to the documentation for that.
 
 ### Fit the model
 
-We now have everything together to specify and fit our model. First, we specify our optimizer backend as
+Now we can construct the *SemOptimizer* that will use the *NLopt* engine for constrained optimization.
 
 ```@example constraints
 using NLopt
 
-constrained_optimizer = SemOptimizerNLopt(
+constrained_optimizer = SemOptimizer(
+    engine = :NLopt,
     algorithm = :AUGLAG,
     options = Dict(:upper_bounds => upper_bounds, :xtol_abs => 1e-4),
     local_algorithm = :LD_LBFGS,
-    equality_constraints = NLoptConstraint(;f = eq_constraint, tol = 1e-8),
-    inequality_constraints = NLoptConstraint(;f = ineq_constraint, tol = 1e-8),
+    equality_constraints = (f => 1e-8),
+    inequality_constraints = (g => 1e-8),
 )
 ```
 
-As you see, the equality constraints and inequality constraints are passed as keyword arguments, and the bounds are passed as options for the (outer) optimization algorithm.
+As you see, the equality and inequality constraints are passed as keyword arguments, and the bounds are passed as options for the (outer) optimization algorithm.
 Additionally, for equality and inequality constraints, a feasibility tolerance can be specified that controls if a solution can be accepted, even if it violates the constraints by a small amount.
 Especially for equality constraints, it is recommended to allow for a small positive tolerance.
 In this example, we set both tolerances to `1e-8`.
 
 !!! warning "Convergence criteria"
     We have often observed that the default convergence criteria in NLopt lead to non-convergence flags.
     Indeed, this example does not convergence with default criteria.
-    As you see above, we used a realively liberal absolute tolerance in the optimization parameters of 1e-4.
+    As you see above, we used a relatively liberal absolute tolerance in the optimization parameters of 1e-4.
     This should not be a problem in most cases, as the sampling variance in (almost all) structural equation models
     should lead to uncertainty in the parameter estimates that are orders of magnitude larger.
     We nontheless recommend choosing a convergence criterion with care (i.e. w.r.t. the scale of your parameters),
     inspecting the solutions for plausibility, and comparing them to unconstrained solutions.
 
-```@example constraints
-model_constrained = Sem(
-    specification = partable,
-    data = data
-)
+We now have everything to fit our model under constraints:
 
-model_fit_constrained = fit(constrained_optimizer, model_constrained)
+```@example constraints
+model_fit_constrained = fit(constrained_optimizer, model)
 ```
 
 As you can see, the optimizer converged (`:XTOL_REACHED`) and investigating the solution yields

diff --git a/docs/src/tutorials/construction/build_by_parts.md b/docs/src/tutorials/construction/build_by_parts.md
@@ -59,7 +59,7 @@ ml = SemML(observed = observed)
 loss_ml = SemLoss(ml)
 
 # optimizer ----------------------------------------------------------------------------
-optimizer = SemOptimizerOptim()
+optimizer = SemOptimizer()
 
 # model --------------------------------------------------------------------------------
 

diff --git a/docs/src/tutorials/construction/outer_constructor.md b/docs/src/tutorials/construction/outer_constructor.md
@@ -41,7 +41,7 @@ model = Sem(
     data = data,
     implied = RAMSymbolic,
     loss = SemWLS,
-    optimizer = SemOptimizerOptim
+    optimizer = SemOptimizer
 )
 ```
 

diff --git a/docs/src/tutorials/fitting/fitting.md b/docs/src/tutorials/fitting/fitting.md
@@ -17,7 +17,6 @@ Structural Equation Model
 - Fields
    observed:  SemObservedData
    implied:   RAM
-   optimizer: SemOptimizerOptim
 
 ------------- Optimization result -------------
 
@@ -60,7 +59,7 @@ The available keyword arguments are listed in the sections [Using Optim.jl](@ref
 Alternative, you can also explicitely define a `SemOptimizer` and pass it as the first argument to `fit`:
 
 ```julia
-my_optimizer = SemOptimizerOptim(algorithm = BFGS())
+my_optimizer = SemOptimizer(algorithm = BFGS())
 
 fit(my_optimizer, model)
 ```