Exact penalty method

exact_penalty_method(M, f, grad_f, p=rand(M); kwargs...)
exact_penalty_method(M, cmo::ConstrainedManifoldObjective, p=rand(M); kwargs...)
exact_penalty_method!(M, f, grad_f, p; kwargs...)
exact_penalty_method!(M, cmo::ConstrainedManifoldObjective, p; kwargs...)

perform the exact penalty method (EPM) [LB19] The aim of the EPM is to find a solution of the constrained optimization task

\[\begin{aligned} \operatorname*{arg\,min}_{p ∈ \mathcal{M}} & f(p)\\ \text{subject to}\quad&g_i(p) ≤ 0 \quad \text{ for } i= 1, …, m,\\ \quad & h_j(p)=0 \quad \text{ for } j=1,…,n, \end{aligned}\]

where M is a Riemannian manifold, and $f$, $\{g_i\}_{i=1}^{m}$ and $\{h_j\}_{j=1}^{n}$ are twice continuously differentiable functions from M to ℝ. For that a weighted $L_1$-penalty term for the violation of the constraints is added to the objective

\[f(p) + ρ\biggl( \sum_{i=1}^m \max\bigl\{0, g_i(p)\bigr\} + \sum_{j=1}^n \vert h_j(p)\vert\biggr),\]

where $ρ>0$ is the penalty parameter.

Since this is non-smooth, a SmoothingTechnique with parameter u is applied, see the ExactPenaltyCost.

In every step $k$ of the exact penalty method, the smoothed objective is then minimized over all $p ∈\mathcal{M}$. Then, the accuracy tolerance $ϵ$ and the smoothing parameter $u$ are updated by setting

\[ϵ^{(k)}=\max\{ϵ_{\min}, θ_ϵ ϵ^{(k-1)}\},\]

where $ϵ_{\min}$ is the lowest value $ϵ$ is allowed to become and $θ_ϵ ∈ (0,1)$ is a constant scaling factor, and

\[u^{(k)} = \max \{u_{\min}, \theta_u u^{(k-1)} \},\]

where $u_{\min}$ is the lowest value $u$ is allowed to become and $θ_u ∈ (0,1)$ is a constant scaling factor.

Finally, the penalty parameter $ρ$ is updated as

\[ρ^{(k)} = \begin{cases} ρ^{(k-1)}/θ_ρ, & \text{if } \displaystyle \max_{j ∈ \mathcal{E},i ∈ \mathcal{I}} \Bigl\{ \vert h_j(p^{(k)}) \vert, g_i(p^{(k)})\Bigr\} > u^{(k-1)},\\ ρ^{(k-1)}, & \text{ else,} \end{cases}\]

where $θ_ρ ∈ (0,1)$ is a constant scaling factor.

Input

• M::AbstractManifold: a Riemannian manifold $\mathcal{M}$
• f: a cost function $f: \mathcal{M}→ ℝ$ implemented as (M, p) -> v
• grad_f: the (Riemannian) gradient $\operatorname{grad}f: \mathcal{M} → T\mathcal{M}$ of f as a function (M, p) -> X or a function (M, X, p) -> X computing X in-place
• p::P: a point on the manifold $\mathcal{M}$

Keyword arguments

if not called with the ConstrainedManifoldObjective cmo

• g=missing: the inequality constraints
• grad_g=missing: the gradient of the inequality constraints
• grad_h=missing: the gradient of the equality constraints
• h=missing: the equality constraints

Note that one of the pairs (g, grad_g) or (h, grad_h) has to be provided. Otherwise the problem is not constrained and a better solver would be for example quasi_Newton.

Further keyword arguments

• callbacks::D = Dict{Symbol,Function}(): provided callback functions given either as a single function (symbol, problem, state, k) called in every hook or as a (vector of) pairs :hook => function, which are processed by process_callbacks_arg. As key you can either pass single symbol or an array of symbols to indicate a callback should be added in multiple places
• equality_constraints=nothing: the number $n$ of equality constraints. If not provided, a call to the gradient of h is performed to estimate these.
• gradient_equality_range=gradient_range: specify how gradients of the equality constraints are represented, see VectorGradientFunction.
• gradient_inequality_range=gradient_range: specify how gradients of the inequality constraints are represented, see VectorGradientFunction.
• gradient_range=nothing: specify how both gradients of the constraints are represented
• inequality_constraints=nothing: the number $m$ of inequality constraints. If not provided, a call to the gradient of g is performed to estimate these.
• smoothing=LogarithmicSumOfExponentials: a SmoothingTechnique to use
• stopping_criterion::StoppingCriterion=StopAfterIteration(300)|(StopWhenSmallerOrEqual(:ϵ, ϵ_min)&StopWhenChangeLess(1e-10) ): a functor indicating that the stopping criterion is fulfilled
• sub_cost=ExactPenaltyCost(cmo, ρ, u; smoothing=smoothing): cost to use in the sub solver. This is used to define the sub_problem= keyword and has hence no effect if you set sub_problem directly.
• sub_grad=ExactPenaltyGrad(cmo, ρ, u; smoothing=smoothing): gradient to use in the sub solver. This is used to define the sub_problem= keyword and has hence no effect if you set sub_problem directly.
• sub_kwargs = (;): a named tuple of keyword arguments that are passed to decorate_objective! of the sub solver's objective, the decorate_state! of the sub solver's state, and the sub state constructor itself.
• sub_problem::Union{AbstractManoptProblem, F} =DefaultManoptProblem(M,ManifoldGradientObjective(sub_cost, sub_grad; evaluation=evaluation)): specify a problem for a solver or a closed form solution function, which can be allocating or in-place.
• sub_state::Union{AbstractManoptSolverState,AbstractEvaluationType} =QuasiNewtonState: a state to specify the sub solver to use. For a closed form solution, this indicates the type of function. The default uses a QuasiNewtonLimitedMemoryDirectionUpdate with InverseBFGS.
• sub_stopping_criterion=StopAfterIteration(300)|StopWhenGradientNormLess(ϵ)|StopWhenStepsizeLess(1e-8): a stopping criterion for the sub solver This is used to define the sub_state= keyword and has hence no effect if you set sub_state directly.
• u=1e-1: the smoothing parameter and threshold for violation of the constraints
• u_exponent=1/100: exponent of the u update factor;
• u_min=1e-6: the lower bound for the smoothing parameter and threshold for violation of the constraints
• ρ=1.0: the penalty parameter
• ϵ=1e-3: the accuracy tolerance
• ϵ_exponent=1/100: exponent of the ϵ update factor;
• ϵ_min=1e-6: the lower bound for the accuracy tolerance

For the ranges of the constraints' gradient, other power manifold tangent space representations, mainly the ArrayPowerRepresentation can be used if the gradients can be computed more efficiently in that representation.

All other keyword arguments are passed to decorate_state! for state decorators or decorate_objective! for objective decorators, respectively.

Output

The obtained approximate minimizer $p^*$. To obtain the whole final state of the solver, see get_solver_return for details, especially the return_state= keyword.

exact_penalty_method(M, f, grad_f, p=rand(M); kwargs...)
exact_penalty_method(M, cmo::ConstrainedManifoldObjective, p=rand(M); kwargs...)
exact_penalty_method!(M, f, grad_f, p; kwargs...)
exact_penalty_method!(M, cmo::ConstrainedManifoldObjective, p; kwargs...)

perform the exact penalty method (EPM) [LB19] The aim of the EPM is to find a solution of the constrained optimization task

\[\begin{aligned} \operatorname*{arg\,min}_{p ∈ \mathcal{M}} & f(p)\\ \text{subject to}\quad&g_i(p) ≤ 0 \quad \text{ for } i= 1, …, m,\\ \quad & h_j(p)=0 \quad \text{ for } j=1,…,n, \end{aligned}\]

where M is a Riemannian manifold, and $f$, $\{g_i\}_{i=1}^{m}$ and $\{h_j\}_{j=1}^{n}$ are twice continuously differentiable functions from M to ℝ. For that a weighted $L_1$-penalty term for the violation of the constraints is added to the objective

\[f(p) + ρ\biggl( \sum_{i=1}^m \max\bigl\{0, g_i(p)\bigr\} + \sum_{j=1}^n \vert h_j(p)\vert\biggr),\]

where $ρ>0$ is the penalty parameter.

Since this is non-smooth, a SmoothingTechnique with parameter u is applied, see the ExactPenaltyCost.

In every step $k$ of the exact penalty method, the smoothed objective is then minimized over all $p ∈\mathcal{M}$. Then, the accuracy tolerance $ϵ$ and the smoothing parameter $u$ are updated by setting

\[ϵ^{(k)}=\max\{ϵ_{\min}, θ_ϵ ϵ^{(k-1)}\},\]

where $ϵ_{\min}$ is the lowest value $ϵ$ is allowed to become and $θ_ϵ ∈ (0,1)$ is a constant scaling factor, and

\[u^{(k)} = \max \{u_{\min}, \theta_u u^{(k-1)} \},\]

where $u_{\min}$ is the lowest value $u$ is allowed to become and $θ_u ∈ (0,1)$ is a constant scaling factor.

Finally, the penalty parameter $ρ$ is updated as

\[ρ^{(k)} = \begin{cases} ρ^{(k-1)}/θ_ρ, & \text{if } \displaystyle \max_{j ∈ \mathcal{E},i ∈ \mathcal{I}} \Bigl\{ \vert h_j(p^{(k)}) \vert, g_i(p^{(k)})\Bigr\} > u^{(k-1)},\\ ρ^{(k-1)}, & \text{ else,} \end{cases}\]

where $θ_ρ ∈ (0,1)$ is a constant scaling factor.

Input

• M::AbstractManifold: a Riemannian manifold $\mathcal{M}$
• f: a cost function $f: \mathcal{M}→ ℝ$ implemented as (M, p) -> v
• grad_f: the (Riemannian) gradient $\operatorname{grad}f: \mathcal{M} → T\mathcal{M}$ of f as a function (M, p) -> X or a function (M, X, p) -> X computing X in-place
• p::P: a point on the manifold $\mathcal{M}$

Keyword arguments

if not called with the ConstrainedManifoldObjective cmo

• g=missing: the inequality constraints
• grad_g=missing: the gradient of the inequality constraints
• grad_h=missing: the gradient of the equality constraints
• h=missing: the equality constraints

Note that one of the pairs (g, grad_g) or (h, grad_h) has to be provided. Otherwise the problem is not constrained and a better solver would be for example quasi_Newton.

Further keyword arguments

• callbacks::D = Dict{Symbol,Function}(): provided callback functions given either as a single function (symbol, problem, state, k) called in every hook or as a (vector of) pairs :hook => function, which are processed by process_callbacks_arg. As key you can either pass single symbol or an array of symbols to indicate a callback should be added in multiple places
• equality_constraints=nothing: the number $n$ of equality constraints. If not provided, a call to the gradient of h is performed to estimate these.
• gradient_equality_range=gradient_range: specify how gradients of the equality constraints are represented, see VectorGradientFunction.
• gradient_inequality_range=gradient_range: specify how gradients of the inequality constraints are represented, see VectorGradientFunction.
• gradient_range=nothing: specify how both gradients of the constraints are represented
• inequality_constraints=nothing: the number $m$ of inequality constraints. If not provided, a call to the gradient of g is performed to estimate these.
• smoothing=LogarithmicSumOfExponentials: a SmoothingTechnique to use
• stopping_criterion::StoppingCriterion=StopAfterIteration(300)|(StopWhenSmallerOrEqual(:ϵ, ϵ_min)&StopWhenChangeLess(1e-10) ): a functor indicating that the stopping criterion is fulfilled
• sub_cost=ExactPenaltyCost(cmo, ρ, u; smoothing=smoothing): cost to use in the sub solver. This is used to define the sub_problem= keyword and has hence no effect if you set sub_problem directly.
• sub_grad=ExactPenaltyGrad(cmo, ρ, u; smoothing=smoothing): gradient to use in the sub solver. This is used to define the sub_problem= keyword and has hence no effect if you set sub_problem directly.
• sub_kwargs = (;): a named tuple of keyword arguments that are passed to decorate_objective! of the sub solver's objective, the decorate_state! of the sub solver's state, and the sub state constructor itself.
• sub_problem::Union{AbstractManoptProblem, F} =DefaultManoptProblem(M,ManifoldGradientObjective(sub_cost, sub_grad; evaluation=evaluation)): specify a problem for a solver or a closed form solution function, which can be allocating or in-place.
• sub_state::Union{AbstractManoptSolverState,AbstractEvaluationType} =QuasiNewtonState: a state to specify the sub solver to use. For a closed form solution, this indicates the type of function. The default uses a QuasiNewtonLimitedMemoryDirectionUpdate with InverseBFGS.
• sub_stopping_criterion=StopAfterIteration(300)|StopWhenGradientNormLess(ϵ)|StopWhenStepsizeLess(1e-8): a stopping criterion for the sub solver This is used to define the sub_state= keyword and has hence no effect if you set sub_state directly.
• u=1e-1: the smoothing parameter and threshold for violation of the constraints
• u_exponent=1/100: exponent of the u update factor;
• u_min=1e-6: the lower bound for the smoothing parameter and threshold for violation of the constraints
• ρ=1.0: the penalty parameter
• ϵ=1e-3: the accuracy tolerance
• ϵ_exponent=1/100: exponent of the ϵ update factor;
• ϵ_min=1e-6: the lower bound for the accuracy tolerance

For the ranges of the constraints' gradient, other power manifold tangent space representations, mainly the ArrayPowerRepresentation can be used if the gradients can be computed more efficiently in that representation.

All other keyword arguments are passed to decorate_state! for state decorators or decorate_objective! for objective decorators, respectively.

Output

The obtained approximate minimizer $p^*$. To obtain the whole final state of the solver, see get_solver_return for details, especially the return_state= keyword.

ExactPenaltyMethodState{P,T} <: AbstractManoptSolverState

Describes the exact penalty method, with

Fields

• callbacks::D: provided callback functions given as a dictionary with symbols as keys
• ϵ: the accuracy tolerance
• ϵ_min: the lower bound for the accuracy tolerance
• p::P: a point on the manifold $\mathcal{M}$ storing the current iterate
• ρ: the penalty parameter
• sub_problem::Union{AbstractManoptProblem, F}: specify a problem for a solver or a closed form solution function, which can be allocating or in-place.
• sub_state::Union{AbstractManoptSolverState,AbstractEvaluationType}: a state to specify the sub solver to use. For a closed form solution, this indicates the type of function.
• stop::StoppingCriterion: a functor indicating that the stopping criterion is fulfilled
• u: the smoothing parameter and threshold for violation of the constraints
• u_min: the lower bound for the smoothing parameter and threshold for violation of the constraints
• θ_ϵ: the scaling factor of the tolerance parameter
• θ_ρ: the scaling factor of the penalty parameter
• θ_u: the scaling factor of the smoothing parameter

Constructor

ExactPenaltyMethodState(M::AbstractManifold, sub_problem, sub_state; kwargs...)

construct the exact penalty state.

ExactPenaltyMethodState(M::AbstractManifold, sub_problem;    evaluation=AllocatingEvaluation(), kwargs...)

construct the exact penalty state, where sub_problem is a closed form solution with evaluation as type of evaluation. The closed form solution is expected to be of the form (M, q, ρ, u, p) -> q for the in-place and (M, ρ, u, p) -> q for the allocating evaluation, that is it minimizes the smoothed penalized function for the current penalty parameter ρ and smoothing parameter u, starting from p.

Keyword arguments

• callbacks::D = Dict{Symbol,Function}(): provided callback functions given as a dictionary with symbols as keys
• p::P =rand(M): a point on the manifold $\mathcal{M}$ to specify the initial value
• u=1e-1
• u_exponent=1 / 100: a shortcut for the scaling factor $θ_u$.
• u_min=1e-6
• θ_u=(u_min / u)^(u_exponent)
• θ_ρ=0.3
• ρ=1.0
• ϵ=1e-3
• ϵ_exponent=1 / 100: a shortcut for the scaling factor $θ_ϵ$
• ϵ_min=1e-6
• stopping_criterion::StoppingCriterion=StopAfterIteration(300)|(StopWhenSmallerOrEqual(:ϵ, ϵ_min)&StopWhenChangeLess(1e-10) ): a functor indicating that the stopping criterion is fulfilled
• θ_ϵ=(ϵ_min / ϵ)^(ϵ_exponent)

See also

exact_penalty_method