diff --git a/Project.toml b/Project.toml
index 37634708..0cc596b9 100644
--- 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
index 1bb68c4d..cf273c05 100644
--- 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
index 9e01ac87..1a094a88 100644
--- 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) = ...
 ```
\ No newline at end of file
diff --git a/docs/src/performance/simulation.md b/docs/src/performance/simulation.md
index d268853f..85a0c0a0 100644
--- 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
index feb5c8f4..b9dac45e 100644
--- a/docs/src/tutorials/backends/nlopt.md
+++ b/docs/src/tutorials/backends/nlopt.md
@@ -1,31 +1,21 @@
 # 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,
@@ -33,15 +23,14 @@ my_optimizer = SemOptimizerNLopt(;
 )
 ```
 
-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
index cf287e77..00fcbe94 100644
--- 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
index 2b453925..8c2c8850 100644
--- 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]
 ```
\ No newline at end of file
diff --git a/docs/src/tutorials/constraints/constraints.md b/docs/src/tutorials/constraints/constraints.md
index c433240a..32bb6a52 100644
--- 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,29 +117,26 @@ 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`.
@@ -141,19 +144,16 @@ 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
index 680e2880..52e12f30 100644
--- 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
index e2772430..e0c69ef3 100644
--- 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
index d7353c9f..1af03ce8 100644
--- 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)
 ```
diff --git a/docs/src/tutorials/regularization/regularization.md b/docs/src/tutorials/regularization/regularization.md
index 2b2c6df3..17add030 100644
--- a/docs/src/tutorials/regularization/regularization.md
+++ b/docs/src/tutorials/regularization/regularization.md
@@ -5,7 +5,9 @@
 For ridge regularization, you can simply use `SemRidge` as an additional loss function
 (for example, a model with the loss functions `SemML` and `SemRidge` corresponds to ridge-regularized maximum likelihood estimation).
 
-For lasso, elastic net and (far) beyond, you can load the `ProximalAlgorithms.jl` and `ProximalOperators.jl` packages alongside `StructuralEquationModels`:
+You can define lasso, elastic net and other forms of regularization using [`ProximalOperators.jl`](https://github.com/JuliaFirstOrder/ProximalOperators.jl)
+and optimize the SEM model with [`ProximalAlgorithms.jl`](https://github.com/JuliaFirstOrder/ProximalAlgorithms.jl)
+that provides so-called *proximal optimization* algorithms.
 
 ```@setup reg
 using StructuralEquationModels, ProximalAlgorithms, ProximalOperators
@@ -19,24 +21,22 @@ Pkg.add("ProximalOperators")
 using StructuralEquationModels, ProximalAlgorithms, ProximalOperators
 ```
 
-## `SemOptimizerProximal`
+## Proximal optimization
 
-To estimate regularized models, we provide a "building block" for the optimizer part, called `SemOptimizerProximal`.
-It connects our package to the [`ProximalAlgorithms.jl`](https://github.com/JuliaFirstOrder/ProximalAlgorithms.jl) optimization backend, providing so-called proximal optimization algorithms.
-Those can handle, amongst other things, various forms of regularization.
-
-It can be used as
+With *ProximalAlgorithms* package loaded, it is now possible to use `:Proximal` optimization engine
+in `SemOptimizer` for estimating regularized models.
 
 ```julia
-SemOptimizerProximal(
+SemOptimizer(;
+    engine = :Proximal,
     algorithm = ProximalAlgorithms.PANOC(),
     operator_g,
     operator_h = nothing
 )
 ```
 
-The proximal operator (aka the regularization function) can be passed as `operator_g`.
-The available Algorithms are listed [here](https://juliafirstorder.github.io/ProximalAlgorithms.jl/stable/guide/implemented_algorithms/).
+The *proximal operator* (aka the *regularization function*) is passed as `operator_g`, see [available operators](https://juliafirstorder.github.io/ProximalOperators.jl/stable/functions/).
+The `algorithm` is chosen from one of the [available algorithms](https://juliafirstorder.github.io/ProximalAlgorithms.jl/stable/guide/implemented_algorithms/).
 
 ## First example - lasso
 
@@ -84,7 +84,7 @@ model = Sem(
 We labeled the covariances between the items because we want to regularize those:
 
 ```@example reg
-ind = getindex.(
+cov_inds = getindex.(
     Ref(param_indices(model)),
     [:cov_15, :cov_24, :cov_26, :cov_37, :cov_48, :cov_68])
 ```
@@ -96,30 +96,24 @@ The lasso penalty is defined as
 \sum \lambda_i \lvert \theta_i \rvert
 ```
 
-From the previously linked [documentation](https://juliafirstorder.github.io/ProximalOperators.jl/stable/functions/#ProximalOperators.NormL1), we find that lasso regularization is named `NormL1` in the `ProximalOperators` package, and that we can pass an array of hyperparameters (`λ`) to control the amount of regularization for each parameter. To regularize only the observed item covariances, we define `λ` as
+In `ProximalOperators.jl`, lasso regularization is represented by the [`NormL1`](https://juliafirstorder.github.io/ProximalOperators.jl/stable/functions/#ProximalOperators.NormL1) operator. It allows controlling the amount of
+regularization individually for each SEM model parameter via the vector of hyperparameters (`λ`).
+To regularize only the observed item covariances, we define `λ` as
 
 ```@example reg
-λ = zeros(31); λ[ind] .= 0.02
-```
-
-and use `SemOptimizerProximal`.
+λ = zeros(31); λ[cov_inds] .= 0.02
 
-```@example reg
-optimizer_lasso = SemOptimizerProximal(
+optimizer_lasso = SemOptimizer(
+    engine = :Proximal,
     operator_g = NormL1(λ)
 )
-
-model_lasso = Sem(
-    specification = partable,
-    data = data
-)
 ```
 
 Let's fit the regularized model
 
 ```@example reg
 
-fit_lasso = fit(optimizer_lasso, model_lasso)
+fit_lasso = fit(optimizer_lasso, model)
 ```
 
 and compare the solution to unregularizted estimates:
@@ -134,34 +128,31 @@ update_partable!(partable, :estimate_lasso, fit_lasso, solution(fit_lasso))
 details(partable)
 ```
 
-Instead of explicitely defining a `SemOptimizerProximal` object, you can also pass `engine = :Proximal` and additional keyword arguments to `fit`:
+Instead of explicitly defining a `SemOptimizer` object, you can also pass `engine = :Proximal`
+and additional keyword arguments directly to the `fit` function:
 
 ```@example reg
-sem_fit = fit(model; engine = :Proximal, operator_g = NormL1(λ))
+fit_lasso2 = fit(model; engine = :Proximal, operator_g = NormL1(λ))
 ```
 
 ## Second example - mixed l1 and l0 regularization
 
 You can choose to penalize different parameters with different types of regularization functions.
-Let's use the lasso again on the covariances, but additionally penalyze the error variances of the observed items via l0 regularization.
+Let's use the *lasso* (*l1*) again on the covariances, but additionally penalize the error variances of the observed items via *l0* regularization.
 
-The l0 penalty is defined as
+The *l0* penalty is defined as
 ```math
-\lambda \mathrm{nnz}(\theta)
+l_0 = \lambda \mathrm{nnz}(\theta)
 ```
 
-To define a sup of separable proximal operators (i.e. no parameter is penalized twice),
-we can use [`SlicedSeparableSum`](https://juliafirstorder.github.io/ProximalOperators.jl/stable/calculus/#ProximalOperators.SlicedSeparableSum) from the `ProximalOperators` package:
+Since we apply *l1* and *l0* to the disjoint sets of parameters, this regularization could be represented as
+as sum of *separable proximal operators* (i.e. no parameter is penalized twice)
+implemented by the [`SlicedSeparableSum`](https://juliafirstorder.github.io/ProximalOperators.jl/stable/calculus/#ProximalOperators.SlicedSeparableSum) operator:
 
 ```@example reg
-prox_operator = SlicedSeparableSum((NormL0(20.0), NormL1(0.02), NormL0(0.0)), ([ind], [9:11], [vcat(1:8, 12:25)]))
-
-model_mixed = Sem(
-    specification = partable,
-    data = data,
-)
+l0_and_l1_reg = SlicedSeparableSum((NormL0(20.0), NormL1(0.02), NormL0(0.0)), ([cov_inds], [9:11], [vcat(1:8, 12:25)]))
 
-fit_mixed = fit(model_mixed; engine = :Proximal, operator_g = prox_operator)
+fit_mixed = fit(model; engine = :Proximal, operator_g = l0_and_l1_reg)
 ```
 
 Let's again compare the different results:
diff --git a/ext/SEMNLOptExt/NLopt.jl b/ext/SEMNLOptExt/NLopt.jl
index 27bc3003..e8409014 100644
--- a/ext/SEMNLOptExt/NLopt.jl
+++ b/ext/SEMNLOptExt/NLopt.jl
@@ -1,37 +1,115 @@
-Base.convert(
-    ::Type{NLoptConstraint},
-    tuple::NamedTuple{(:f, :tol), Tuple{F, T}},
-) where {F, T} = NLoptConstraint(tuple.f, tuple.tol)
+############################################################################################
+### Types
+############################################################################################
+
+const NLoptConstraint = Pair{Any, Number}
+
+struct SemOptimizerNLopt <: SemOptimizer{:NLopt}
+    algorithm::Symbol
+    local_algorithm::Union{Symbol, Nothing}
+    options::Dict{Symbol, Any}
+    local_options::Dict{Symbol, Any}
+    equality_constraints::Vector{NLoptConstraint}
+    inequality_constraints::Vector{NLoptConstraint}
+end
 
 ############################################################################################
 ### Constructor
 ############################################################################################
 
+"""
+Uses *NLopt.jl* as the optimization engine. For more information on the available algorithms 
+and options, see the [*NLopt.jl*](https://github.com/JuliaOpt/NLopt.jl)  package and the 
+[NLopt docs](https://nlopt.readthedocs.io/en/latest/).
+
+# Constructor
+
+    SemOptimizer(;
+        engine = :NLopt,
+        algorithm = :LD_LBFGS,
+        options = Dict{Symbol, Any}(),
+        local_algorithm = nothing,
+        local_options = Dict{Symbol, Any}(),
+        equality_constraints = nothing,
+        inequality_constraints = nothing,
+        constraint_tol::Number = 0.0,
+        kwargs...)
+
+# Arguments
+- `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: optional equality constraints
+- `inequality_constraints:: optional inequality constraints
+- `constraint_tol::Number`: default tolerance for constraints
+
+## Constraints specification
+
+Equality and inequality constraints arguments could be a single constraint or any
+iterable constraints container (e.g. vector or tuple).
+Each constraint could be a function or any other callable object that
+takes the two input arguments:
+  - the vector of the model parameters;
+  - the array for the in-place calculation of the constraint gradient.
+To override the default tolerance, the constraint can be specified
+as a pair of the function and its tolerance: `constraint_func => tol`.
+For information on how to use inequality and equality constraints,
+see [Constrained optimization](@ref) in our online documentation.
+
+# Example
+```julia
+my_optimizer = SemOptimizer(engine = :NLopt)
+
+# constrained optimization with augmented lagrangian
+my_constrained_optimizer = SemOptimizer(;
+    engine = :NLopt,
+    algorithm = :AUGLAG,
+    local_algorithm = :LD_LBFGS,
+    local_options = Dict(:ftol_rel => 1e-6),
+    inequality_constraints = (my_constraint => tol),
+)
+```
+
+# Interfaces
+- `algorithm(::SemOptimizerNLopt)`
+- `local_algorithm(::SemOptimizerNLopt)`
+- `options(::SemOptimizerNLopt)`
+- `local_options(::SemOptimizerNLopt)`
+- `equality_constraints(::SemOptimizerNLopt)`
+- `inequality_constraints(::SemOptimizerNLopt)`
+"""
 function SemOptimizerNLopt(;
     algorithm = :LD_LBFGS,
     local_algorithm = nothing,
     options = Dict{Symbol, Any}(),
     local_options = Dict{Symbol, Any}(),
-    equality_constraints = Vector{NLoptConstraint}(),
-    inequality_constraints = Vector{NLoptConstraint}(),
-    kwargs...,
+    equality_constraints = nothing,
+    inequality_constraints = nothing,
+    constraint_tol::Number = 0.0,
+    kwargs..., # FIXME remove the sink for unused kwargs
 )
-    applicable(iterate, equality_constraints) && !isa(equality_constraints, NamedTuple) ||
-        (equality_constraints = [equality_constraints])
-    applicable(iterate, inequality_constraints) &&
-    !isa(inequality_constraints, NamedTuple) ||
-        (inequality_constraints = [inequality_constraints])
+    constraint(f::Any) = f => constraint_tol
+    constraint(f_and_tol::Pair) = f_and_tol
+
+    constraints(::Nothing) = Vector{NLoptConstraint}()
+    constraints(constraints) =
+        applicable(iterate, constraints) && !isa(constraints, Pair) ?
+        [constraint(constr) for constr in constraints] : [constraint(constraints)]
+
     return SemOptimizerNLopt(
         algorithm,
         local_algorithm,
         options,
         local_options,
-        convert.(NLoptConstraint, equality_constraints),
-        convert.(NLoptConstraint, inequality_constraints),
+        constraints(equality_constraints),
+        constraints(inequality_constraints),
     )
 end
 
-SEM.SemOptimizer{:NLopt}(args...; kwargs...) = SemOptimizerNLopt(args...; kwargs...)
+SEM.SemOptimizer(::Val{:NLopt}, args...; kwargs...) = SemOptimizerNLopt(args...; kwargs...)
+
+SEM.optimizer_engine_doc(engine::Val{:NLopt}) = doc(SemOptimizerNLopt)
 
 ############################################################################################
 ### Recommended methods
@@ -44,29 +122,29 @@ SEM.update_observed(optimizer::SemOptimizerNLopt, observed::SemObserved; kwargs.
 ### additional methods
 ############################################################################################
 
-SEM.algorithm(optimizer::SemOptimizerNLopt) = optimizer.algorithm
 local_algorithm(optimizer::SemOptimizerNLopt) = optimizer.local_algorithm
 SEM.options(optimizer::SemOptimizerNLopt) = optimizer.options
 local_options(optimizer::SemOptimizerNLopt) = optimizer.local_options
 equality_constraints(optimizer::SemOptimizerNLopt) = optimizer.equality_constraints
 inequality_constraints(optimizer::SemOptimizerNLopt) = optimizer.inequality_constraints
 
-mutable struct NLoptResult
+struct NLoptResult
     result::Any
     problem::Any
 end
 
-SEM.optimizer(res::NLoptResult) = res.problem.algorithm
+SEM.algorithm_name(res::NLoptResult) = res.problem.algorithm
 SEM.n_iterations(res::NLoptResult) = res.problem.numevals
 SEM.convergence(res::NLoptResult) = res.result[3]
 
 # construct SemFit from fitted NLopt object
-function SemFit_NLopt(optimization_result, model::AbstractSem, start_val, opt)
+function SemFit_NLopt(optimization_result, model::AbstractSem, start_val, optim, opt)
     return SemFit(
         optimization_result[1],
         optimization_result[2],
         start_val,
         model,
+        optim,
         NLoptResult(optimization_result, opt),
     )
 end
@@ -78,10 +156,7 @@ function SEM.fit(
     start_params::AbstractVector;
     kwargs...,
 )
-
-    # construct the NLopt problem
-    opt = construct_NLopt_problem(optim.algorithm, optim.options, length(start_params))
-    set_NLopt_constraints!(opt, optim)
+    opt = construct_NLopt(optim.algorithm, optim.options, nparams(model))
     opt.min_objective =
         (par, G) -> SEM.evaluate!(
             zero(eltype(par)),
@@ -90,27 +165,30 @@ function SEM.fit(
             model,
             par,
         )
+    for (f, tol) in optim.inequality_constraints
+        inequality_constraint!(opt, f, tol)
+    end
+    for (f, tol) in optim.equality_constraints
+        equality_constraint!(opt, f, tol)
+    end
 
     if !isnothing(optim.local_algorithm)
-        opt_local = construct_NLopt_problem(
-            optim.local_algorithm,
-            optim.local_options,
-            length(start_params),
-        )
+        opt_local =
+            construct_NLopt(optim.local_algorithm, optim.local_options, nparams(model))
         opt.local_optimizer = opt_local
     end
 
     # fit
     result = NLopt.optimize(opt, start_params)
 
-    return SemFit_NLopt(result, model, start_params, opt)
+    return SemFit_NLopt(result, model, start_params, optim, opt)
 end
 
 ############################################################################################
 ### additional functions
 ############################################################################################
 
-function construct_NLopt_problem(algorithm, options, npar)
+function construct_NLopt(algorithm, options, npar)
     opt = Opt(algorithm, npar)
 
     for (key, val) in pairs(options)
@@ -120,15 +198,6 @@ function construct_NLopt_problem(algorithm, options, npar)
     return opt
 end
 
-function set_NLopt_constraints!(opt::Opt, optimizer::SemOptimizerNLopt)
-    for con in optimizer.inequality_constraints
-        inequality_constraint!(opt, con.f, con.tol)
-    end
-    for con in optimizer.equality_constraints
-        equality_constraint!(opt, con.f, con.tol)
-    end
-end
-
 ############################################################################################
 # pretty printing
 ############################################################################################
diff --git a/ext/SEMNLOptExt/SEMNLOptExt.jl b/ext/SEMNLOptExt/SEMNLOptExt.jl
index bf905e3a..b0091dc7 100644
--- a/ext/SEMNLOptExt/SEMNLOptExt.jl
+++ b/ext/SEMNLOptExt/SEMNLOptExt.jl
@@ -1,10 +1,13 @@
 module SEMNLOptExt
 
 using StructuralEquationModels, NLopt
-import StructuralEquationModels: SemOptimizerNLopt, NLoptConstraint
+
+import Base.Docs: doc
 
 SEM = StructuralEquationModels
 
+export SemOptimizerNLopt
+
 include("NLopt.jl")
 
 end
diff --git a/ext/SEMProximalOptExt/ProximalAlgorithms.jl b/ext/SEMProximalOptExt/ProximalAlgorithms.jl
index 0d4748e3..b9bc18bb 100644
--- a/ext/SEMProximalOptExt/ProximalAlgorithms.jl
+++ b/ext/SEMProximalOptExt/ProximalAlgorithms.jl
@@ -1,5 +1,34 @@
+############################################################################################
+### Types
+############################################################################################
+mutable struct SemOptimizerProximal{A, B, C} <: SemOptimizer{:Proximal}
+    algorithm::A
+    operator_g::B
+    operator_h::C
+end
+
 SEM.SemOptimizer{:Proximal}(args...; kwargs...) = SemOptimizerProximal(args...; kwargs...)
 
+SEM.optimizer_engine_doc(engine::Val{:Proximal}) = doc(SemOptimizerProximal)
+
+"""
+Connects to `ProximalAlgorithms.jl` as the optimization backend. For more information on 
+the available algorithms and options, see the online docs on [Regularization](@ref) and
+the documentation of [*ProximalAlgorithms.jl*](https://github.com/JuliaFirstOrder/ProximalAlgorithms.jl) / [ProximalOperators.jl](https://github.com/JuliaFirstOrder/ProximalOperators.jl).
+
+# Constructor
+
+    SemOptimizerProximal(;
+        algorithm = ProximalAlgorithms.PANOC(),
+        operator_g,
+        operator_h = nothing,
+        kwargs...,
+
+# Arguments
+- `algorithm`: optimization algorithm.
+- `operator_g`: proximal operator (e.g., regularization penalty)
+- `operator_h`: optional second proximal operator
+"""
 SemOptimizerProximal(;
     algorithm = ProximalAlgorithms.PANOC(),
     operator_g,
@@ -14,19 +43,12 @@ SemOptimizerProximal(;
 SEM.update_observed(optimizer::SemOptimizerProximal, observed::SemObserved; kwargs...) =
     optimizer
 
-############################################################################################
-### additional methods
-############################################################################################
-
-SEM.algorithm(optimizer::SemOptimizerProximal) = optimizer.algorithm
-
 ############################################################################
-### Pretty Printing
+### Model fitting
 ############################################################################
 
-function Base.show(io::IO, struct_inst::SemOptimizerProximal)
-    print_type_name(io, struct_inst)
-    print_field_types(io, struct_inst)
+mutable struct ProximalResult
+    result::Any
 end
 
 ## connect to ProximalAlgorithms.jl
@@ -36,10 +58,6 @@ function ProximalAlgorithms.value_and_gradient(model::AbstractSem, params)
     return obj, grad
 end
 
-mutable struct ProximalResult
-    result::Any
-end
-
 function SEM.fit(
     optim::SemOptimizerProximal,
     model::AbstractSem,
@@ -75,14 +93,31 @@ function SEM.fit(
         solution,
         start_params,
         model,
+        optim,
         ProximalResult(optimization_result),
     )
 end
 
+############################################################################################
+### additional methods
+############################################################################################
+
+SEM.algorithm_name(res::ProximalResult) = SEM.algorithm_name(res.result[:algorithm])
+SEM.algorithm_name(::ProximalAlgorithms.IterativeAlgorithm{I,H,S,D,K}) where 
+    {I, H, S, D, K} = nameof(I)
+
+SEM.convergence(::ProximalResult) = "No standard convergence criteria for proximal \n algorithms available."
+SEM.n_iterations(res::ProximalResult) = res.result[:iterations]
+
 ############################################################################################
 # pretty printing
 ############################################################################################
 
+function Base.show(io::IO, struct_inst::SemOptimizerProximal)
+    print_type_name(io, struct_inst)
+    print_field_types(io, struct_inst)
+end
+
 function Base.show(io::IO, result::ProximalResult)
     print(io, "Minimum:          $(round(result.result[:minimum]; digits = 2)) \n")
     print(io, "No. evaluations:  $(result.result[:iterations]) \n")
diff --git a/ext/SEMProximalOptExt/SEMProximalOptExt.jl b/ext/SEMProximalOptExt/SEMProximalOptExt.jl
index 04be35cb..2b35fb0c 100644
--- a/ext/SEMProximalOptExt/SEMProximalOptExt.jl
+++ b/ext/SEMProximalOptExt/SEMProximalOptExt.jl
@@ -1,9 +1,11 @@
 module SEMProximalOptExt
 
-using StructuralEquationModels
+using StructuralEquationModels, ProximalAlgorithms
 using StructuralEquationModels: print_type_name, print_field_types
-using ProximalAlgorithms
-import StructuralEquationModels: SemOptimizerProximal
+
+import Base.Docs: doc
+
+export SemOptimizerProximal
 
 SEM = StructuralEquationModels
 
diff --git a/src/StructuralEquationModels.jl b/src/StructuralEquationModels.jl
index 46692bd5..d3372f27 100644
--- a/src/StructuralEquationModels.jl
+++ b/src/StructuralEquationModels.jl
@@ -18,6 +18,10 @@ using LinearAlgebra,
 
 import StatsAPI: params, coef, coefnames, dof, fit, nobs, coeftable
 
+import Base.Docs: doc
+
+using InteractiveUtils: subtypes
+
 export StenoGraphs, @StenoGraph, meld, SimpleNode
 
 const SEM = StructuralEquationModels
@@ -86,9 +90,6 @@ include("frontend/fit/fitmeasures/fit_measures.jl")
 # standard errors
 include("frontend/fit/standard_errors/hessian.jl")
 include("frontend/fit/standard_errors/bootstrap.jl")
-# extensions
-include("package_extensions/SEMNLOptExt.jl")
-include("package_extensions/SEMProximalOptExt.jl")
 
 export AbstractSem,
     AbstractSemSingle,
@@ -123,8 +124,8 @@ export AbstractSem,
     SemWLS,
     loss,
     SemOptimizer,
-    SemOptimizerEmpty,
-    SemOptimizerOptim,
+    optimizer_engine_doc,
+    optimizer_engines,
     optimizer,
     n_iterations,
     convergence,
@@ -196,8 +197,5 @@ export AbstractSem,
     →,
     ←,
     ↔,
-    ⇔,
-    SemOptimizerNLopt,
-    NLoptConstraint,
-    SemOptimizerProximal
+    ⇔
 end
diff --git a/src/frontend/fit/SemFit.jl b/src/frontend/fit/SemFit.jl
index 438da4da..9ed73653 100644
--- a/src/frontend/fit/SemFit.jl
+++ b/src/frontend/fit/SemFit.jl
@@ -13,15 +13,16 @@ Fitted structural equation model.
 - `model(::SemFit)`
 - `optimization_result(::SemFit)`
 
-- `optimizer(::SemFit)` -> optimization algorithm
+- `algorithm_name(::SemFit)` -> optimization algorithm
 - `n_iterations(::SemFit)` -> number of iterations
 - `convergence(::SemFit)` -> convergence properties
 """
-mutable struct SemFit{Mi, So, St, Mo, O}
+mutable struct SemFit{Mi, So, St, Mo, Op, O}
     minimum::Mi
     solution::So
     start_val::St
     model::Mo
+    optimizer::Op
     optimization_result::O
 end
 
@@ -39,6 +40,10 @@ function Base.show(io::IO, semfit::SemFit)
     #print(io, "Objective value: $(round(semfit.minimum, digits = 4)) \n")
     print(io, "------------- Optimization result ------------- \n")
     print(io, "\n")
+    print(io, "engine: ")
+    print(io, optimizer_engine(semfit))
+    print(io, "\n")
+    print(io, "\n")
     print(io, semfit.optimization_result)
 end
 
@@ -58,6 +63,7 @@ model(sem_fit::SemFit) = sem_fit.model
 optimization_result(sem_fit::SemFit) = sem_fit.optimization_result
 
 # optimizer properties
-optimizer(sem_fit::SemFit) = optimizer(optimization_result(sem_fit))
+algorithm_name(sem_fit::SemFit) = algorithm_name(sem_fit.optimization_result)
+optimizer_engine(sem_fit::SemFit) = optimizer_engine(sem_fit.optimizer)
 n_iterations(sem_fit::SemFit) = n_iterations(optimization_result(sem_fit))
 convergence(sem_fit::SemFit) = convergence(optimization_result(sem_fit))
diff --git a/src/frontend/fit/summary.jl b/src/frontend/fit/summary.jl
index 3071d565..435b1747 100644
--- a/src/frontend/fit/summary.jl
+++ b/src/frontend/fit/summary.jl
@@ -7,7 +7,8 @@ function details(sem_fit::SemFit; show_fitmeasures = false, color = :light_cyan,
         color = color,
     )
     print("\n")
-    println("Optimization algorithm:      $(optimizer(sem_fit))")
+    println("Optimization engine:         $(optimizer_engine(sem_fit))")
+    println("Optimization algorithm:      $(algorithm_name(sem_fit))")
     println("Convergence:                 $(convergence(sem_fit))")
     println("No. iterations/evaluations:  $(n_iterations(sem_fit))")
     print("\n")
diff --git a/src/optimizer/Empty.jl b/src/optimizer/Empty.jl
index 1bf0c30a..e7d027df 100644
--- a/src/optimizer/Empty.jl
+++ b/src/optimizer/Empty.jl
@@ -1,13 +1,10 @@
 ############################################################################################
 ### Types
 ############################################################################################
-"""
-Empty placeholder for models that don't need
-an optimizer part.
-
-# Constructor
 
-    SemOptimizerEmpty()
+# dummy SEM optimizer
+"""
+Test.
 """
 struct SemOptimizerEmpty <: SemOptimizer{:Empty} end
 
@@ -15,7 +12,9 @@ struct SemOptimizerEmpty <: SemOptimizer{:Empty} end
 ### Constructor
 ############################################################################################
 
-SemOptimizer{:Empty}() = SemOptimizerEmpty()
+SemOptimizer(::Val{:Empty}) = SemOptimizerEmpty()
+
+optimizer_engine_doc(engine::Val{:Empty}) = doc(SemOptimizerEmpty)
 
 ############################################################################################
 ### Recommended methods
diff --git a/src/optimizer/abstract.jl b/src/optimizer/abstract.jl
index c1ad7259..b05d8c20 100644
--- a/src/optimizer/abstract.jl
+++ b/src/optimizer/abstract.jl
@@ -1,3 +1,79 @@
+const optimizer_engine_packages = Dict(
+    :NLopt => "NLopt",
+    :Proximal => "ProximalAlgorithms"
+)
+
+function throw_engine_error(E)
+    if typeof(E) !== Symbol
+        throw(ArgumentError("engine argument must be a Symbol."))
+    elseif haskey(optimizer_engine_packages, E)
+        error("optimizer \":$E\" requires \"using $(optimizer_engine_packages[E])\".")
+    else
+        error("optimizer engine \":$E\" is not supported.")
+    end
+end
+
+"""
+    SemOptimizer(args...; engine::Symbol = :Optim, kwargs...)
+
+Constructs a `SemOptimizer` object that can be passed to [`fit`](@ref) for specifying aspects 
+of the numerical optimization involved in fitting a SEM.
+
+The keyword `engine` controlls which Julia package is used, with `:Optim` being the default.
+- `optimizer_engines()` prints a list of currently available engines.
+- `optimizer_engine_doc(EngineName)` prints information on the usage of a specific engine.
+
+More engines become available if specific packages are loaded, for example 
+[*NLopt.jl*](https://github.com/JuliaOpt/NLopt.jl) (also see [Constrained optimization](@ref) 
+in the online documentation) or 
+[*ProximalAlgorithms.jl*](https://github.com/JuliaFirstOrder/ProximalAlgorithms.jl)
+(also see [Regularization](@ref) in the online documentation).
+
+The additional arguments `args...` and `kwargs...` are engine-specific and control further
+aspects of the optimization process, such as the algorithm, convergence criteria or constraints.
+Information on those can be accessed with `optimizer_engine_doc`.
+
+To connect the SEM package to a completely new optimization backend, you can implement a new 
+subtype of SemOptimizer.
+"""
+SemOptimizer(args...; engine::Symbol = :Optim, kwargs...) =
+    SemOptimizer{engine}(args...; kwargs...)
+
+# fallback optimizer constructor when the engine E is not supported
+SemOptimizer(::Val{E}, args...; kwargs...) where {E} = throw_engine_error(E)
+
+SemOptimizer{E}(args...; kwargs...) where {E} = SemOptimizer(Val(E), args...; kwargs...)
+
+"""
+    (1) optimizer_engine(::Type{<:SemOptimizer{E}})
+    (2) optimizer_engine(::SemOptimizer{E})
+    
+Returns `E`; the engine of a `SemOptimizer` object or a subtype of `SemOptimizer`.
+"""
+optimizer_engine(::Type{<:SemOptimizer{E}}) where {E} = E
+optimizer_engine(::SemOptimizer{E}) where {E} = E
+
+"""
+    optimizer_engines()
+
+Returns a vector of optimizer engines supported by the `engine` keyword argument of
+the [`SemOptimizer`](@ref) constructor.
+
+The list of engines depends on the Julia packages loaded (with the `using` directive)
+into the current session.
+"""
+optimizer_engines() = Symbol[optimizer_engine(opt_type) for opt_type in subtypes(SemOptimizer)]
+
+"""
+    optimizer_engine_doc(engine::Symbol)
+
+Shows information on the optimizer engine.
+For a list of available engines, call `optimizer_engines`.
+"""
+optimizer_engine_doc(engine) = optimizer_engine_doc(Val(engine))
+
+optimizer_engine_doc(::Val{E}) where {E} = throw_engine_error(E)
+
 """
     fit([optim::SemOptimizer], model::AbstractSem;
             [engine::Symbol], start_val = start_val, kwargs...)
diff --git a/src/optimizer/optim.jl b/src/optimizer/optim.jl
index 2d782473..6a5f617e 100644
--- a/src/optimizer/optim.jl
+++ b/src/optimizer/optim.jl
@@ -3,34 +3,36 @@
 ############################################################################################
 ### Types and Constructor
 ############################################################################################
-"""
-    SemOptimizerOptim{A, B} <: SemOptimizer{:Optim}
 
-Connects to `Optim.jl` as the optimization backend.
+# SemOptimizer for the Optim.jl
+mutable struct SemOptimizerOptim{A, B} <: SemOptimizer{:Optim}
+    algorithm::A
+    options::B
+end
+
+"""
+Connects to *Optim.jl* as the optimization engine.
+For more information on the available algorithms and options, see the [*Optim.jl* docs](https://julianlsolvers.github.io/Optim.jl/stable/).
 
 # Constructor
 
-    SemOptimizerOptim(;
+    SemOptimizer(;
+        engine = :Optim,
         algorithm = LBFGS(),
         options = Optim.Options(;f_reltol = 1e-10, x_abstol = 1.5e-8),
         kwargs...)
 
 # Arguments
-- `algorithm`: optimization algorithm from `Optim.jl`
+- `algorithm`: optimization algorithm from *Optim.jl*
 - `options::Optim.Options`: options for the optimization algorithm
 
-# Usage
-All algorithms and options from the Optim.jl library are available, for more information see
-the Optim.jl online documentation.
-
 # Examples
 ```julia
-my_optimizer = SemOptimizerOptim()
-
 # hessian based optimization with backtracking linesearch and modified initial step size
 using Optim, LineSearches
 
-my_newton_optimizer = SemOptimizerOptim(
+my_newton_optimizer = SemOptimizer(
+    engine = :Optim,
     algorithm = Newton(
         ;linesearch = BackTracking(order=3),
         alphaguess = InitialHagerZhang()
@@ -38,10 +40,7 @@ my_newton_optimizer = SemOptimizerOptim(
 )
 ```
 
-# Extended help
-
-## Constrained optimization
-
+# Constrained optimization
 When using the `Fminbox` or `SAMIN` constrained optimization algorithms,
 the vector or dictionary of lower and upper bounds for each model parameter can be specified
 via `lower_bounds` and `upper_bounds` keyword arguments.
@@ -50,27 +49,20 @@ the default bound for all non-variance model parameters,
 and the `variance_lower_bound` and `variance_upper_bound` keyword --
 for the variance parameters (the diagonal of the *S* matrix).
 
-## Interfaces
-- `algorithm(::SemOptimizerOptim)`
-- `options(::SemOptimizerOptim)`
-
-## Implementation
-
-Subtype of `SemOptimizer`.
+# Interfaces
+- `algorithm(::SemOptimizer)`
+- `options(::SemOptimizer)`
 """
-mutable struct SemOptimizerOptim{A, B} <: SemOptimizer{:Optim}
-    algorithm::A
-    options::B
-end
-
-SemOptimizer{:Optim}(args...; kwargs...) = SemOptimizerOptim(args...; kwargs...)
-
 SemOptimizerOptim(;
     algorithm = LBFGS(),
     options = Optim.Options(; f_reltol = 1e-10, x_abstol = 1.5e-8),
     kwargs...,
 ) = SemOptimizerOptim(algorithm, options)
 
+SemOptimizer(::Val{:Optim}, args...; kwargs...) = SemOptimizerOptim(args...; kwargs...)
+
+SEM.optimizer_engine_doc(engine::Val{:Optim}) = doc(SemOptimizerOptim)
+
 ############################################################################################
 ### Recommended methods
 ############################################################################################
@@ -81,24 +73,9 @@ update_observed(optimizer::SemOptimizerOptim, observed::SemObserved; kwargs...)
 ### additional methods
 ############################################################################################
 
-algorithm(optimizer::SemOptimizerOptim) = optimizer.algorithm
 options(optimizer::SemOptimizerOptim) = optimizer.options
 
-function SemFit(
-    optimization_result::Optim.MultivariateOptimizationResults,
-    model::AbstractSem,
-    start_val,
-)
-    return SemFit(
-        optimization_result.minimum,
-        optimization_result.minimizer,
-        start_val,
-        model,
-        optimization_result,
-    )
-end
-
-optimizer(res::Optim.MultivariateOptimizationResults) = Optim.summary(res)
+algorithm_name(res::Optim.MultivariateOptimizationResults) = Optim.summary(res)
 n_iterations(res::Optim.MultivariateOptimizationResults) = Optim.iterations(res)
 convergence(res::Optim.MultivariateOptimizationResults) = Optim.converged(res)
 
@@ -147,5 +124,11 @@ function fit(
             optim.options,
         )
     end
-    return SemFit(result, model, start_params)
+    return SemFit(
+        result.minimum,
+        result.minimizer,
+        start_params,
+        model,
+        optim,
+        result)
 end
diff --git a/src/package_extensions/SEMNLOptExt.jl b/src/package_extensions/SEMNLOptExt.jl
deleted file mode 100644
index 64c4cff0..00000000
--- a/src/package_extensions/SEMNLOptExt.jl
+++ /dev/null
@@ -1,65 +0,0 @@
-"""
-Connects to `NLopt.jl` as the optimization backend.
-Only usable if `NLopt.jl` is loaded in the current Julia session!
-
-# Constructor
-
-    SemOptimizerNLopt(;
-        algorithm = :LD_LBFGS,
-        options = Dict{Symbol, Any}(),
-        local_algorithm = nothing,
-        local_options = Dict{Symbol, Any}(),
-        equality_constraints = Vector{NLoptConstraint}(),
-        inequality_constraints = Vector{NLoptConstraint}(),
-        kwargs...)
-
-# Arguments
-- `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
-
-# Example
-```julia
-my_optimizer = SemOptimizerNLopt()
-
-# constrained optimization with augmented lagrangian
-my_constrained_optimizer = SemOptimizerNLopt(;
-    algorithm = :AUGLAG,
-    local_algorithm = :LD_LBFGS,
-    local_options = Dict(:ftol_rel => 1e-6),
-    inequality_constraints = NLoptConstraint(;f = my_constraint, tol = 0.0),
-)
-```
-
-# Usage
-All algorithms and options from the NLopt library are available, for more information see
-the NLopt.jl package and the NLopt online documentation.
-For information on how to use inequality and equality constraints,
-see [Constrained optimization](@ref) in our online documentation.
-
-# Extended help
-
-## Interfaces
-- `algorithm(::SemOptimizerNLopt)`
-- `local_algorithm(::SemOptimizerNLopt)`
-- `options(::SemOptimizerNLopt)`
-- `local_options(::SemOptimizerNLopt)`
-- `equality_constraints(::SemOptimizerNLopt)`
-- `inequality_constraints(::SemOptimizerNLopt)`
-"""
-struct SemOptimizerNLopt{A, A2, B, B2, C} <: SemOptimizer{:NLopt}
-    algorithm::A
-    local_algorithm::A2
-    options::B
-    local_options::B2
-    equality_constraints::C
-    inequality_constraints::C
-end
-
-Base.@kwdef struct NLoptConstraint
-    f::Any
-    tol = 0.0
-end
diff --git a/src/package_extensions/SEMProximalOptExt.jl b/src/package_extensions/SEMProximalOptExt.jl
deleted file mode 100644
index ad4c2da2..00000000
--- a/src/package_extensions/SEMProximalOptExt.jl
+++ /dev/null
@@ -1,27 +0,0 @@
-"""
-Connects to `ProximalAlgorithms.jl` as the optimization backend.
-Can be used for regularized SEM, for a tutorial see the online docs on [Regularization](@ref).
-
-# Constructor
-
-    SemOptimizerProximal(;
-        algorithm = ProximalAlgorithms.PANOC(),
-        operator_g,
-        operator_h = nothing,
-        kwargs...,
-
-# Arguments
-- `algorithm`: optimization algorithm.
-- `operator_g`: proximal operator (e.g., regularization penalty)
-- `operator_h`: optional second proximal operator
-
-# Usage
-All algorithms and operators from `ProximalAlgorithms.jl` are available,
-for more information see the online docs on [Regularization](@ref) and
-the documentation of `ProximalAlgorithms.jl` / `ProximalOperators.jl`.
-"""
-mutable struct SemOptimizerProximal{A, B, C} <: SemOptimizer{:Proximal}
-    algorithm::A
-    operator_g::B
-    operator_h::C
-end
diff --git a/src/types.jl b/src/types.jl
index 0e279e5b..92ca2c28 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -79,24 +79,8 @@ end
 Base.:*(x::SemWeight{Nothing}, y) = y
 Base.:*(x::SemWeight, y) = x.w * y
 
-"""
-Supertype of all objects that can serve as the `optimizer` field of a SEM.
-Connects the SEM to its optimization backend and controls options like the optimization algorithm.
-If you want to connect the SEM package to a new optimization backend, you should implement a subtype of SemOptimizer.
-"""
 abstract type SemOptimizer{E} end
 
-engine(::Type{SemOptimizer{E}}) where {E} = E
-engine(optimizer::SemOptimizer) = engine(typeof(optimizer))
-
-SemOptimizer(args...; engine::Symbol = :Optim, kwargs...) =
-    SemOptimizer{engine}(args...; kwargs...)
-
-# fallback optimizer constructor
-function SemOptimizer{E}(args...; kwargs...) where {E}
-    throw(ErrorException("$E optimizer is not supported."))
-end
-
 """
 Supertype of all objects that can serve as the observed field of a SEM.
 Pre-processes data and computes sufficient statistics for example.
diff --git a/test/examples/multigroup/multigroup.jl b/test/examples/multigroup/multigroup.jl
index 43de554c..2d43c3d2 100644
--- a/test/examples/multigroup/multigroup.jl
+++ b/test/examples/multigroup/multigroup.jl
@@ -86,7 +86,7 @@ start_test = [
     fill(0.05, 3)
     fill(0.01, 3)
 ]
-semoptimizer = SemOptimizerOptim()
+semoptimizer = SemOptimizer()
 
 @testset "RAMMatrices | constructor | Optim" begin
     include("build_models.jl")
@@ -169,7 +169,7 @@ start_test = [
     0.01
     0.05
 ]
-semoptimizer = SemOptimizerOptim()
+semoptimizer = SemOptimizer()
 
 @testset "Graph → Partable → RAMMatrices | constructor | Optim" begin
     include("build_models.jl")
diff --git a/test/examples/political_democracy/constraints.jl b/test/examples/political_democracy/constraints.jl
index cc1b0874..7a6670fa 100644
--- a/test/examples/political_democracy/constraints.jl
+++ b/test/examples/political_democracy/constraints.jl
@@ -26,8 +26,8 @@ constrained_optimizer = SemOptimizer(;
     algorithm = :AUGLAG,
     local_algorithm = :LD_LBFGS,
     options = Dict(:xtol_rel => 1e-4),
-    # equality_constraints = (f = eq_constraint, tol = 1e-14),
-    inequality_constraints = (f = ineq_constraint, tol = 0.0),
+    # equality_constraints = (eq_constraint => 1e-14),
+    inequality_constraints = (ineq_constraint => 0.0),
 )
 
 @test constrained_optimizer isa SemOptimizer{:NLopt}
diff --git a/test/examples/recover_parameters/recover_parameters_twofact.jl b/test/examples/recover_parameters/recover_parameters_twofact.jl
index ce7dc61f..9f9503af 100644
--- a/test/examples/recover_parameters/recover_parameters_twofact.jl
+++ b/test/examples/recover_parameters/recover_parameters_twofact.jl
@@ -68,7 +68,7 @@ loss_ml = SemLoss(SemML(; observed = semobserved, nparams = length(start)))
 model_ml = Sem(semobserved, implied_ml, loss_ml)
 objective!(model_ml, true_val)
 
-optimizer = SemOptimizerOptim(
+optimizer = SemOptimizer(
     BFGS(; linesearch = BackTracking(order = 3), alphaguess = InitialHagerZhang()),# m = 100),
     Optim.Options(; f_reltol = 1e-10, x_abstol = 1.5e-8),
 )