Interior point Newton method

interior_point_Newton(M, f, grad_f, Hess_f, p=rand(M); kwargs...)
interior_point_Newton(M, cmo::ConstrainedManifoldObjective, p=rand(M); kwargs...)
interior_point_Newton!(M, f, grad]_f, Hess_f, p; kwargs...)
interior_point_Newton!(M, cmo::ConstrainedManifoldObjective, p; kwargs...)

perform the interior point Newton method following [LY24].

In order to solve the constrained problem

\[\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}\]

This algorithms iteratively solves the linear system based on extending the KKT system by a slack variable s.

\[ \operatorname{J} F(p, μ, λ, s)[X, Y, Z, W] = -F(p, μ, λ, s), \text{ where } X ∈ T_{p}\mathcal{M}, Y,W ∈ ℝ^m, Z ∈ ℝ^n\]

see CondensedKKTVectorFieldJacobian and CondensedKKTVectorField, respectively, for the reduced form, this is usually solved in. From the resulting X and Z in the reduced form, the other two, $Y$, $W$, are then computed.

From the gradient $(X,Y,Z,W)$ at the current iterate $(p, μ, λ, s)$, a line search is performed using the KKTVectorFieldNormSq norm of the KKT vector field (squared) and its gradient KKTVectorFieldNormSqGradient together with the InteriorPointCentralityCondition.

Note that since the vector field $F$ includes the gradients of the constraint functions $g, h$, its gradient or Jacobian requires the Hessians of the constraints.

For that search direction a line search is performed, that additionally ensures that the constraints are further fulfilled.

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_{p}\mathcal{M}$ of f as a function (M, p) -> X or a function (M, X, p) -> X computing X in-place
• Hess_f: the (Riemannian) Hessian $\operatorname{Hess}f: T_{p}\mathcal{M} → T_{p}\mathcal{M}$ of f as a function (M, p, X) -> Y or a function (M, Y, p, X) -> Y computing Y in-place
• p::P: a point on the manifold $\mathcal{M}$

or a ConstrainedManifoldObjective cmo containing f, grad_f, Hess_f, and the constraints

Keyword arguments

The keyword arguments related to the constraints (the first eleven) are ignored if you pass a ConstrainedManifoldObjective cmo

• callbacks::D = Dict{Symbol,Function}(): provided callback functions . A single function (symbol, problem, state, k) called in every hook, a (vector of) pairs :hook => function or will be 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
• centrality_condition=missing; an additional condition when to accept a step size. This can be used to ensure that the resulting iterate is still an interior point if you provide a check (N,q) -> true/false, where N is the manifold of the step_problem.
• equality_constraints=nothing: the number $n$ of equality constraints.
• evaluation::AbstractEvaluationType=AllocatingEvaluation(): specify whether the functions that return an array, for example a point or a tangent vector, work by allocating its result (AllocatingEvaluation) or whether they modify their input argument to return the result therein (InplaceEvaluation). Since usually the first argument is the manifold, the modified argument is the second.
• g=missing: the inequality constraints
• grad_g=missing: the gradient of the inequality constraints
• grad_h=missing: the gradient of the equality constraints
• gradient_range=nothing: specify how gradients are represented, where nothing is equivalent to NestedPowerRepresentation
• gradient_equality_range=gradient_range: specify how the gradients of the equality constraints are represented
• gradient_inequality_range=gradient_range: specify how the gradients of the inequality constraints are represented
• h=missing: the equality constraints
• Hess_g=missing: the Hessian of the inequality constraints
• Hess_h=missing: the Hessian of the equality constraints
• inequality_constraints=nothing: the number $m$ of inequality constraints.
• λ=ones(length(h(M, p))): the Lagrange multiplier with respect to the equality constraints $h$
• μ=ones(length(g(M, p))): the Lagrange multiplier with respect to the inequality constraints $g$
• retraction_method::AbstractRetractionMethod=default_retraction_method(M, typeof(p)): a retraction $\operatorname{retr}$ to use, see the section on retractions
• ρ=μ's / length(μ): store the orthogonality μ's/m to compute the barrier parameter β in the sub problem.
• s=copy(μ): initial value for the slack variables
• σ=calculate_σ(M, cmo, p, μ, λ, s): scaling factor for the barrier parameter β in the sub problem, which is updated during the iterations
• step_objective: a ManifoldGradientObjective of the norm of the KKT vector field KKTVectorFieldNormSq and its gradient KKTVectorFieldNormSqGradient
• step_problem: the manifold $\mathcal{M} × ℝ^m × ℝ^n × ℝ^m$ together with the step_objective as the problem the line search stepsize= employs for determining a step size
• step_state: the StepsizeState with point and search direction
• stepsize::Stepsize=ArmijoLinesearch(): a functor inheriting from Stepsize to determine a step size with the centrality_condition keyword as additional criterion to accept a step, if this is provided"))
• stopping_criterion::StoppingCriterion=StopAfterIteration(200)|StopWhenKKTResidualLess(1e-8): a functor indicating that the stopping criterion is fulfilled a stopping criterion, by default depending on the residual of the KKT vector field or a maximal number of steps, which ever hits first.
• sub_kwargs=(;): keyword arguments to decorate the sub options, for example debug, that automatically respects the main solvers debug options (like sub-sampling) as well
• sub_objective: The SymmetricLinearSystemObjective modelling the system of equations to use in the sub solver, includes the CondensedKKTVectorFieldJacobian $\mathcal{A}(X)$ and the CondensedKKTVectorField $b$ in $\mathcal{A}(X) + b = 0$ we aim to solve. This is used to define the sub_problem= keyword and has hence no effect, if you set sub_problem directly.
• sub_stopping_criterion=StopAfterIteration(manifold_dimension(M))|StopWhenRelativeResidualLess(c,1e-8), where $c = \lVert b \rVert$ from the system to solve. This is used to define the sub_state= keyword and has hence no effect, if you set sub_state directly.
• sub_problem::Union{AbstractManoptProblem, F} =DefaultManoptProblem(M, sub_objective): specify a problem for a solver or a closed form solution function, which can be allocating or in-place.
• sub_state::Union{AbstractManoptProblem, F} =ConjugateResidualState: a state to specify the sub solver to use. For a closed form solution, this indicates the type of function.
• vector_space=Rn a function that, given an integer, returns the manifold to be used for the vector space components $ℝ^m,ℝ^n$
• X=zero_vector(M,p): the initial gradient with respect to p.
• Y=zero(μ): the initial gradient with respect to μ
• Z=zero(λ): the initial gradient with respect to λ
• W=zero(s): the initial gradient with respect to s
• is_feasible_error=:error: specify how to handle infeasible starting points, see is_feasible for options.

As well as internal keywords used to set up these given keywords like _step_M, _step_p, _sub_M, _sub_p, and _sub_X, that should not be changed.

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

Output

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

interior_point_Newton(M, f, grad_f, Hess_f, p=rand(M); kwargs...)
interior_point_Newton(M, cmo::ConstrainedManifoldObjective, p=rand(M); kwargs...)
interior_point_Newton!(M, f, grad]_f, Hess_f, p; kwargs...)
interior_point_Newton!(M, cmo::ConstrainedManifoldObjective, p; kwargs...)

perform the interior point Newton method following [LY24].

In order to solve the constrained problem

\[\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}\]

This algorithms iteratively solves the linear system based on extending the KKT system by a slack variable s.

\[ \operatorname{J} F(p, μ, λ, s)[X, Y, Z, W] = -F(p, μ, λ, s), \text{ where } X ∈ T_{p}\mathcal{M}, Y,W ∈ ℝ^m, Z ∈ ℝ^n\]

see CondensedKKTVectorFieldJacobian and CondensedKKTVectorField, respectively, for the reduced form, this is usually solved in. From the resulting X and Z in the reduced form, the other two, $Y$, $W$, are then computed.

From the gradient $(X,Y,Z,W)$ at the current iterate $(p, μ, λ, s)$, a line search is performed using the KKTVectorFieldNormSq norm of the KKT vector field (squared) and its gradient KKTVectorFieldNormSqGradient together with the InteriorPointCentralityCondition.

Note that since the vector field $F$ includes the gradients of the constraint functions $g, h$, its gradient or Jacobian requires the Hessians of the constraints.

For that search direction a line search is performed, that additionally ensures that the constraints are further fulfilled.

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_{p}\mathcal{M}$ of f as a function (M, p) -> X or a function (M, X, p) -> X computing X in-place
• Hess_f: the (Riemannian) Hessian $\operatorname{Hess}f: T_{p}\mathcal{M} → T_{p}\mathcal{M}$ of f as a function (M, p, X) -> Y or a function (M, Y, p, X) -> Y computing Y in-place
• p::P: a point on the manifold $\mathcal{M}$

or a ConstrainedManifoldObjective cmo containing f, grad_f, Hess_f, and the constraints

Keyword arguments

The keyword arguments related to the constraints (the first eleven) are ignored if you pass a ConstrainedManifoldObjective cmo

• callbacks::D = Dict{Symbol,Function}(): provided callback functions . A single function (symbol, problem, state, k) called in every hook, a (vector of) pairs :hook => function or will be 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
• centrality_condition=missing; an additional condition when to accept a step size. This can be used to ensure that the resulting iterate is still an interior point if you provide a check (N,q) -> true/false, where N is the manifold of the step_problem.
• equality_constraints=nothing: the number $n$ of equality constraints.
• evaluation::AbstractEvaluationType=AllocatingEvaluation(): specify whether the functions that return an array, for example a point or a tangent vector, work by allocating its result (AllocatingEvaluation) or whether they modify their input argument to return the result therein (InplaceEvaluation). Since usually the first argument is the manifold, the modified argument is the second.
• g=missing: the inequality constraints
• grad_g=missing: the gradient of the inequality constraints
• grad_h=missing: the gradient of the equality constraints
• gradient_range=nothing: specify how gradients are represented, where nothing is equivalent to NestedPowerRepresentation
• gradient_equality_range=gradient_range: specify how the gradients of the equality constraints are represented
• gradient_inequality_range=gradient_range: specify how the gradients of the inequality constraints are represented
• h=missing: the equality constraints
• Hess_g=missing: the Hessian of the inequality constraints
• Hess_h=missing: the Hessian of the equality constraints
• inequality_constraints=nothing: the number $m$ of inequality constraints.
• λ=ones(length(h(M, p))): the Lagrange multiplier with respect to the equality constraints $h$
• μ=ones(length(g(M, p))): the Lagrange multiplier with respect to the inequality constraints $g$
• retraction_method::AbstractRetractionMethod=default_retraction_method(M, typeof(p)): a retraction $\operatorname{retr}$ to use, see the section on retractions
• ρ=μ's / length(μ): store the orthogonality μ's/m to compute the barrier parameter β in the sub problem.
• s=copy(μ): initial value for the slack variables
• σ=calculate_σ(M, cmo, p, μ, λ, s): scaling factor for the barrier parameter β in the sub problem, which is updated during the iterations
• step_objective: a ManifoldGradientObjective of the norm of the KKT vector field KKTVectorFieldNormSq and its gradient KKTVectorFieldNormSqGradient
• step_problem: the manifold $\mathcal{M} × ℝ^m × ℝ^n × ℝ^m$ together with the step_objective as the problem the line search stepsize= employs for determining a step size
• step_state: the StepsizeState with point and search direction
• stepsize::Stepsize=ArmijoLinesearch(): a functor inheriting from Stepsize to determine a step size with the centrality_condition keyword as additional criterion to accept a step, if this is provided"))
• stopping_criterion::StoppingCriterion=StopAfterIteration(200)|StopWhenKKTResidualLess(1e-8): a functor indicating that the stopping criterion is fulfilled a stopping criterion, by default depending on the residual of the KKT vector field or a maximal number of steps, which ever hits first.
• sub_kwargs=(;): keyword arguments to decorate the sub options, for example debug, that automatically respects the main solvers debug options (like sub-sampling) as well
• sub_objective: The SymmetricLinearSystemObjective modelling the system of equations to use in the sub solver, includes the CondensedKKTVectorFieldJacobian $\mathcal{A}(X)$ and the CondensedKKTVectorField $b$ in $\mathcal{A}(X) + b = 0$ we aim to solve. This is used to define the sub_problem= keyword and has hence no effect, if you set sub_problem directly.
• sub_stopping_criterion=StopAfterIteration(manifold_dimension(M))|StopWhenRelativeResidualLess(c,1e-8), where $c = \lVert b \rVert$ from the system to solve. This is used to define the sub_state= keyword and has hence no effect, if you set sub_state directly.
• sub_problem::Union{AbstractManoptProblem, F} =DefaultManoptProblem(M, sub_objective): specify a problem for a solver or a closed form solution function, which can be allocating or in-place.
• sub_state::Union{AbstractManoptProblem, F} =ConjugateResidualState: a state to specify the sub solver to use. For a closed form solution, this indicates the type of function.
• vector_space=Rn a function that, given an integer, returns the manifold to be used for the vector space components $ℝ^m,ℝ^n$
• X=zero_vector(M,p): the initial gradient with respect to p.
• Y=zero(μ): the initial gradient with respect to μ
• Z=zero(λ): the initial gradient with respect to λ
• W=zero(s): the initial gradient with respect to s
• is_feasible_error=:error: specify how to handle infeasible starting points, see is_feasible for options.

As well as internal keywords used to set up these given keywords like _step_M, _step_p, _sub_M, _sub_p, and _sub_X, that should not be changed.

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

Output

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

InteriorPointNewtonState{P,T} <: AbstractHessianSolverState

Fields

• callbacks::D: provided callback functions given as a dictionary with symbols as keys
• λ: the Lagrange multiplier with respect to the equality constraints
• μ: the Lagrange multiplier with respect to the inequality constraints
• p::P: a point on the manifold $\mathcal{M}$ storing the current iterate
• s: the current slack variable
• 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{AbstractManoptProblem, F}: a state to specify the sub solver to use. For a closed form solution, this indicates the type of function.
• X: the current gradient with respect to p
• Y: the current gradient with respect to μ
• Z: the current gradient with respect to λ
• W: the current gradient with respect to s
• ρ: store the orthogonality μ's/m to compute the barrier parameter β in the sub problem
• σ: scaling factor for the barrier parameter β in the sub problem
• stop::StoppingCriterion: a functor indicating that the stopping criterion is fulfilled
• retraction_method::AbstractRetractionMethod: a retraction $\operatorname{retr}$ to use, see the section on retractions
• stepsize::Stepsize: a functor inheriting from Stepsize to determine a step size
• step_problem: an AbstractManoptProblem storing the manifold and objective for the line search
• step_state: storing iterate and search direction in a state for the line search, see StepsizeState

Constructor

InteriorPointNewtonState(    M::AbstractManifold, cmo::ConstrainedManifoldObjective, sub_problem::Pr, sub_state::St; kwargs...)InteriorPointNewtonState(    M::AbstractManifold, cmo::ConstrainedManifoldObjective, sub_problem::Pr; kwargs...)InteriorPointNewtonState(sub_problem::Pr, sub_state::St; kwargs...)

Initialize the state, where both the AbstractManifold and the ConstrainedManifoldObjective are used to fill in reasonable defaults for the keywords. For a closed form solution of the sub solver, you can provide the evaluation either as St in the first constructor or as a keyword like in the second. The third constructor is considered an internal constructor accepting the same keywords, but those that are filled by defaults based on M or cmo become mandatory

Input

• M::AbstractManifold: a Riemannian manifold $\mathcal{M}$
• cmo: a ConstrainedManifoldObjective
• 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{AbstractManoptProblem, F}: a state to specify the sub solver to use. For a closed form solution, this indicates the type of function.

Keyword arguments

Let m and n denote the number of inequality and equality constraints, respectively

• callbacks::D = Dict{Symbol,Function}(): provided callback functions given as a dictionary with symbols as keys
• is_feasible_error=:error: specify how to handle infeasible starting points, see is_feasible for options.
• p::P =rand(M): a point on the manifold $\mathcal{M}$ to specify the initial value
• retraction_method::AbstractRetractionMethod=default_retraction_method(M, typeof(p)): a retraction $\operatorname{retr}$ to use, see the section on retractions
• s=ones(m) slack variables for the inequality constraints
• step_objective=ManifoldGradientObjective(KKTVectorFieldNormSq(cmo), KKTVectorFieldNormSqGradient(cmo); evaluation=InplaceEvaluation())
• step_problem: wrap the manifold $\mathcal{M} × ℝ^m × ℝ^n × ℝ^m$
• step_state: the StepsizeState with point and search direction
• stepsize::Stepsize=ArmijoLinesearch(): a functor inheriting from Stepsize to determine a step size with the InteriorPointCentralityCondition as additional condition to accept a step"))
• stopping_criterion::StoppingCriterion=StopAfterIteration(200)|StopWhenChangeLess(1e-8): a functor indicating that the stopping criterion is fulfilled
• vector_space=Rn: a function that, given an integer, returns the manifold to be used for the vector space components $ℝ^m,ℝ^n$
• W=zero(s) tangent vector (gradient) for the slack variables
• X=zero_vector(M,p)
• Y=zero(μ) tangent vector (gradient) for the inequality constraints
• Z=zero(λ) tangent vector (gradient) for the equality constraints
• λ=zeros(n) Lagrange multipliers for the equality constraints
• μ=ones(m) Lagrange multipliers for the inequality constraints
• ρ=μ's/m storage for the orthogonality check
• σ=calculate_σ(M, cmo, p, μ, λ, s)

and internally _step_M and _step_p for the manifold and point in the stepsize.

StepsizeState{P,T} <: AbstractManoptSolverState

A state to store a point and a descent direction used within a linesearch, if these are different from the iterate and search direction of the main solver.

Fields

• p::P: a point on a manifold
• X::T: a tangent vector at p.

Constructor

StepsizeState(p,X)StepsizeState(M::AbstractManifold; p=rand(M), x=zero_vector(M,p)

See also

interior_point_Newton

InteriorPointCentralityCondition{CO,R}

A functor to check the centrality condition.

In order to obtain a step in the linesearch performed within the interior_point_Newton, Section 6 of [LY24] propose the following additional conditions to hold inspired by the Euclidean case described in Section 6 [ETTZ96]:

For a given ConstrainedManifoldObjective assume consider the KKTVectorField $F$, that is we are at a point $q = (p, λ, μ, s)$ on $\mathcal{M} × ℝ^m × ℝ^n × ℝ^m$and a search direction $V = (X, Y, Z, W)$.

Then, let

\[τ_1 = \frac{m\min\{μ ⊙ s\}}{μ^{\mathrm{T}}s} \quad\text{ and }\quad τ_2 = \frac{μ^{\mathrm{T}}s}{\lVert F(q) \rVert},\]

where $⊙$ denotes the Hadamard (or elementwise) product.

For a new candidate $q(α) = \bigl(p(α), λ(α), μ(α), s(α)\bigr := (\operatorname{retr}_p(αX), λ+αY, μ+αZ, s+αW)$, we then define two functions

\[c_1(α) = \min\{μ(α) ⊙ s(α)\} - \frac{γτ_1 μ(α)^{\mathrm{T}}s(α)}{m} \quad\text{ and }\quad c_2(α) = μ(α)^{\mathrm{T}}s(α) – γτ_2 \lVert F(q(α)) \rVert.\]

While the paper now states that the (Armijo) line search starts at a point $\tilde α$, it is easier to include the condition that $c_1(α) ≥ 0$ and $c_2(α) ≥ 0$ into the line search as well.

The functor InteriorPointCentralityCondition(cmo, γ, μ, s, normKKT)(N,qα) defined here evaluates this condition and returns true if both $c_1$ and $c_2$ are non-negative.

Fields

• cmo: a ConstrainedManifoldObjective
• γ: a constant
• τ1, τ2: the constants given in the formula.

Constructor

InteriorPointCentralityCondition(cmo, γ)InteriorPointCentralityCondition(cmo, γ, τ1, τ2)

Initialize the centrality conditions. The parameters τ1, τ2 are initialize to zero if not provided.

calculate_σ(M, cmo, p, μ, λ, s; kwargs...)

Compute the new $σ$ factor for the barrier parameter in interior_point_Newton as

\[\min\{\frac{1}{2}, \lVert F(p; μ, λ, s) \rVert^{\frac{1}{2}}\},\]

where $F$ is the KKT vector field, hence the KKTVectorFieldNormSq is used.

Keyword arguments

• vector_space=Rn a function that, given an integer, returns the manifold to be used for the vector space components $ℝ^m,ℝ^n$
• N the manifold $\mathcal{M} × ℝ^m × ℝ^n × ℝ^m$ the vector field lives on (generated using vector_space)
• q provide memory on N for interims evaluation of the vector field

StopWhenKKTResidualLess <: StoppingCriterion

Stop when the KKT residual

r^2= \lVert \operatorname{grad}_p \mathcal{L}(p, μ, λ)  \rVert^2+ \sum_{i=1}^{m} [μ_i]_{-}^2 + [g_i(p)]_+^2 + \lvert μ_i g_i(p) \rvert^2+ \sum_{j=1}^{n} \lvert h_i(p) \rvert^2.

is less than a given threshold $r < ε$. We use $[v]_+ = \max\{0,v\}$ and $[v]_- = \min\{0,t\}$ for the positive and negative part of $v$, respectively

Fields

• ε: a threshold
• residual: store the last residual if the stopping criterion is hit.
• at_iteration::Int: an integer indicating at which the stopping criterion last indicted to stop, which might also be before the solver started (0). Any negative value indicates that this was not yet the case;