# Usage

This section provides details on how to instantiate and solve a constrained least-squares problem with Enlsip.jl As a reminder from Home, problems to solve are of the following form:

\begin{aligned} \min_{x \in \mathbb{R}^n} \quad & \dfrac{1}{2} \|r(x)\|^2 \\ \text{s.t.} \quad & c_i(x) = 0, \quad i \in \mathcal{E} \\ & c_i(x) \geq 0, \quad i \in \mathcal{I}, \\ \end{aligned}

with:

• the residuals $r_i:\mathbb{R}^n\rightarrow\mathbb{R}$ and the constraints $c_i:\mathbb{R}^n\rightarrow\mathbb{R}$ assumed to be $\mathcal{C}^2$ functions;
• norm $\|\cdot\|$ denoting the Euclidean norm.

Note that with this formulation, bounds constraints are not distinguished from general inequality constraints. Though, for ease of use, they can be provided directly as vectors of lower and/or upper bounds (see Instantiate).

It should be borne in mind, however, that the method implemented in Enlsip.jl has been conceived for nonlinear problems, as there is no other assumption made about the nature of the residuals and constraints functions, apart from being two-time continuously differentiable. The algorithm can still be used to solve linear least squares subject to linear constraints but it will not be as effective as other software where those aspects are taken into account in the design of the optimization method.

## Instantiate a model

Solving a problem with Enlsip is organized in two steps.

First, a model of type CnlsModel must be instantiated.

The CnlsModel constructor requires the evaluation functions of residuals, constraints, their associated jacobian matrices and dimensions of the problem.

Although the package enables one to create linear unconstrained least squares, it is recommended to use it on problems with nonlinear residuals and general constraints.

The three following positional arguments are mandatory to create a model:

• residuals : function that computes the vector of residuals
• nb_parameters : number of variables
• nb_residuals : number of residuals

The following keyword arguments are optional and deal with constraints and Jacobian matrices. If the Jacobian matrix functions are not provided, they are computed numerically by forward differences using automatic differenciation[Backend].

ArgumentDetails
starting_pointinitial solution (can be an infeasbile point)
jacobian_residualsfunction computing the Jacobian matrix of the residuals
eq_constraintsfunction computing the equality constraints
jacobian_eqconsfunction computing the Jacobian matrix of the equality constraints
nb_eqconsnumber of equality constraints
ineq_constraintsfunction computing the inequality constraints
jacobian_ineqconsfunction computing the Jacobian matrix of the inequality constraints
nb_ineqconsnumber of inequality constraints
x_lowvector of lower bounds
x_uppvector of upper bounds

It is assumed that the the different functions passed as arguments of the CnlsModel constructor are called f(x), where x is a vector of nb_parameters elements and f is one of the functions residuals, eq_constraints, jacobian_eqcons etc.

## Solving

Then, the Enlsip solver can be used by calling the solve! function on an instantiated model. One can pass keywords arguments to modify some options of the algorithm or adjust the tolerances. For the expression of the termination criteria, see Method page (Stopping).

Enlsip.solve!Function
solve!(model{T})

Once a CnlsModel has been instantiated, this function solves the optimzation problem associated by using the method implemented in the Enlsip solver.

Options:

• silent::Bool

• Set to false if one wants the algorithm to print details about the iterations and termination of the solver

• Default is true, i.e. by default, there is no output. If one wants to print those information afert solving, the print_cnls_model method

can be called

• max_iter::Int

• Maximum number of iterations allowed

• Default is 100

• scaling::Bool

• Set to true if one wants the algorithm to work with a constraints jacobian matrix whose rows are scaled (i.e. all constraints gradients vectors are scaled)

• Default is false

• time_limit::T

• Maximum elapsed time (i.e. wall time)

• Default is 1000

Tolerances:

• abs_tol::T

• Absolute tolerance for small residuals

• Default is eps(T)

• rel_tol::T

• Relative tolerance used to measure first order criticality and consistency

• Default is sqrt(eps(T))

• c_tol::T

• Tolerance used to measure feasability of the constraints

• Default is sqrt(eps(T))

• x_tol::T

• Tolerance used to measure the distance between two consecutive iterates

• Default is sqrt(eps(T))

Diagnosis of the conduct of the algorithm can be printed by either setting the silent keyword argument of the function solve! to false or by calling print_cnls_model after solving. Here are some details on how to read and understand the different columns of the output:

ColumnDescription
iteriteration number
objectivevalue of the sum of squared residuals (i.e. objective function) at current point
$\vert\vert$ active_constraints $\vert\vert^2$value of the sum of squared active constraints at current point
$\vert\vert$ p $\vert\vert$norm of the search direction computed at current iteration
$\alpha$value of the steplength computed at current iteration
reductionreduction in the objective function performed after moving to the next iterate

One can get additional info about termination of the algorithm by calling one of the following functions:

Name
solution
status
constraints_values
sum_sq_residuals
Enlsip.solutionFunction
solution(model)

Once the given model has been solved, this function returns the optimal solution, or last solution obtained if no convergence, as a Vector of approriate dimension.

Enlsip.statusFunction
status(model)

This functions returns a Symbol that gives brief information on the solving status of model.

If a model has been instantiated but the solver has not been called yet, it will return :unsolved.

Once the solver has been called and if a first order stationary point satisfying the convergence criteria has been computed, it will return :found_first_order_stationary_point.

If the algorithm met an abnormall termination criteria, it will return one of the following:

• :failed : the algorithm encoutered a numerical error that triggered termination

• :maximum_iterations_exceeded : a solution could not be reached within the maximum number of iterations

• :time_limit_exceeded : the algorithm stopped because solving time exceeded the time limit

Enlsip.constraints_valuesFunction
constraints_values(model)

Computes values of all the constraints in model at the solution.

The vector returned is the concatenation of equalities, inequalities and box constraints (in that order).

For instance, let xₛ be the solution found. If functions h, g compute equality and inequality constraints and xₗ, xᵤ are vectors of lower and lower bounds, it will return [h(xₛ); g(xₛ); xₛ-xₗ; xᵤ-xₛ].

If one wants to compute each type of constraints seperately, see equality_constraints_values, inequality_constraints_values and bounds_constraints_values.

Enlsip.sum_sq_residualsFunction
sum_sq_residuals(model)

Once the given model has been solved, returns the value of the objective function, i.e. sum of squared residuals functions, computed at the optimal solution. If no convergence, this value is computed at the last solution obtained.

## Examples

### Problem 65 from Hock and Schittkowski collection[HS80]

We show how to implement and solve the following problem:

\begin{aligned} \min_{x_1, x_2, x_3} \quad & (x_1-x_2)^2 + \dfrac{(x_1+x_2-10)^2}{9}+(x_3-5)^2 \\ \text{s.t.} \quad & 48-x_1^2-x_2^2-x_3^2 \geq 0\\ & -4.5\leq x_i \leq 4.5, \quad i=1,2\\ & -5 \leq x_3 \leq 5. \end{aligned}.

The expected optimal solution is $(3.650461821, 3.65046168, 4.6204170507)$.

Associated value of objective function equals $0.9535288567$.

First, we provide the dimensions of the problems.

# Dimensions of the problem

n = 3 # number of parameters
m = 3 # number of residuals

Then, we define the functions required to compute the residuals, constraints, their respective jacobian matrices and a starting point. In this example, we use the starting point given in the reference[HS80], i.e. $(-5, 5, 0)$

# Residuals and Jacobian matrix associated
r(x::Vector) = [x[1] - x[2]; (x[1]+x[2]-10.0) / 3.0; x[3]-5.0]

jac_r(x::Vector) = [1. -1. 0;
1/3 1/3 0.;
0. 0. 1.]

# Constraints (one nonlinear inequality and box constraints)

c(x::Vector) = [48.0 - x[1]^2-x[2]^2-x[3]^2] # evaluation function for the equality constraint
jac_c(x::Vector) = [ -2x[1] -2x[2] -2x[3]] # Jacobian matrix of the equality constraint

x_l = [-4.5, -4.5, -5.0] # lower bounds
x_u = [4.5, 4.5, 5.0] # upper bounds

# Starting point
x0 = [-5.0, 5.0, 0.0]

A CnlsModel can now be instantiated.

# Instantiate a model associated with the problem
hs65_model = Enlsip.CnlsModel(r, n, m ;starting_point=x0, ineq_constraints = c,
nb_ineqcons = 1, x_low=x_l, x_upp=x_u, jacobian_residuals=jac_r, jacobian_ineqcons=jac_c)

Finally, the solve! function can be called on our model. In this example, keyword arguments remain to default values.

Enlsip.solve!(hs65_model)

Once Enlsip solver has been executed on our problem, a summary of the conduct of the algorithm can be printed by calling print_cnls_model.

Enlsip.print_cnls_model(hs65_model)


****************************************************************
*                                                              *
*                       Enlsip.jl v0.9.7                       *
*                                                              *
* This is the Julia version of the ENLSIP algorithm, initially *
* conceived and developed in Fortran77 by Per Lindstrom and    *
* Per Ake Wedin from the Institute of Information Processing   *
* (University of Umea, Sweden).                                *
*                                                              *
****************************************************************

Characteristics of the model

Number of parameters.................:     3
Number of residuals..................:     3
Number of equality constraints.......:     0
Number of inequality constraints.....:     1
Number of lower bounds...............:     3
Number of upper bounds...............:     3
Constraints internal scaling.........: false

Iteration steps information

iter    objective   ||active_constraints||²  ||p||       α     reduction
1  1.3611111e+02       4.50e+00         5.00e+00  4.81e-01  1.599e+01
2  1.1589067e+02       4.66e+01         6.84e+00  4.30e-01  2.657e+01
3  7.8319097e+01       6.61e-02         8.19e+00  7.58e-01  7.366e+01
4  4.6650091e+00       2.18e+01         1.44e+00  1.00e+00  5.459e+00
5  9.5479609e-01       4.29e+00         4.66e-01  1.00e+00  4.411e-01
6  9.3767054e-01       4.70e-02         6.67e-02  1.00e+00  1.766e-02
7  9.5320179e-01       1.98e-05         7.92e-03  1.00e+00  3.653e-04
8  9.5352428e-01       3.94e-09         9.97e-04  1.00e+00  5.152e-06
9  9.5352878e-01       9.89e-13         1.25e-04  1.00e+00  8.164e-08
10  9.5352886e-01       2.46e-16         1.57e-05  1.00e+00  1.287e-09
11  9.5352886e-01       6.11e-20         1.98e-06  1.00e+00  2.030e-11
12  9.5352886e-01       1.52e-23         2.48e-07  1.00e+00  3.205e-13
13  9.5352886e-01       3.65e-27         3.12e-08  1.00e+00  4.774e-15

Number of iterations...................:   13

Square sum of residuals................: 9.5352886e-01

Number of function evaluations.........:  152
Number of Jacobian matrix evaluations..:   30

Solving time (seconds).................: 0.080
Termination status.....................: found_first_order_stationary_point

If one just wants to know about termination of the algorithm, calling status will tell if the problem has been successfully solved or not.

Enlsip.status(hs65_model)
:found_first_order_stationary_point

Then, calling solution and sum_sq_residuals will respectively return the optimal solution obtained and the value of objective function at that point.

hs65_solution = Enlsip.solution(hs65_model)
3-element Vector{Float64}:
3.6504617268193496
3.650461726819378
4.620417552781781
hs65_objective = Enlsip.sum_sq_residuals(hs65_model)
0.9535288568047824

The solution obtained is relatively close to the expected optimal solution, although it differs from more than the tolerance used.

maximum(abs.(hs65_solution - [3.650461821, 3.65046168, 4.6204170507])) < sqrt(eps(Float64))
false

However, the difference between the objective value obtained with Enlsip and the expected one does not exceed the default tolerance.

abs(hs65_objective - 0.9535288567) < sqrt(eps(Float64))
true

### Chained Rosenbrock problem[LV99]

We now consider the following problem:

\begin{aligned} \min_{x} \quad & \sum_{i=1}^{n-1} 100(x_i^2-x_{i+1})^2 + (x_i-1)^2 \\ \text{s.t.} \quad & 3x_{k+1}^3 + 2x_{k+2}+\sin(x_{k+1}-x_{k+2})\sin(x_{k+1}+x_{k+2})+4x_{k+1} - x_k\exp(x_k-x_{k+1}) - 8 = 0 \\ & k=1,\ldots, n-2, \end{aligned}

for a given natural number $n\geq 3$. In this example, we use $n=1000$. Though analytic Jacobians are easy to define, they will be passed as arguments in the model instantiation, so to let the algorithm use automatic differentiation.

# Dimensions
n = 1000 # Number of variables
m = 2(n-1) # Number of residuals
nb_eq = n-2 # Number of equality constraints

# Residuals
function r(x::Vector)
n = length(x)
m = 2(n-1)
rx = Vector{eltype(x)}(undef,m)
rx[1:n-1] = [10(x[i]^2 - x[i+1]) for i=1:n-1]
rx[n:m] = [x[k-n+1] - 1 for k=n:m]
return rx
end

# Constraints
function c(x::Vector)
n = length(x)
cx = [3x[k+1]^3 + 2x[k+2] - 5 + sin(x[k+1]-x[k+2])*sin(x[k+1]+x[k+2]) + 4x[k+1] -
x[k]*exp(x[k]-x[k+1]) - 3 for k=1:n-2]
return cx
end

x0 = [(mod(i,2) == 1 ? -1.2 : 1.0) for i=1:n] # Starting point

# Instantiation of the model
cr_enlsip = CnlsModel(r,n,m; starting_point=x0, eq_constraints=c, nb_eqcons=nb_eq)

# Solving
Enlsip.solve!(cr_enlsip)

# Show solving status
Enlsip.status(cr_enlsip)
:found_first_order_stationary_point

To give an insight on how the performance of the algorithm scales when dimensions increase, we give a table of the solving time obtained with different values of parameter n. For comparison purposes, the calculation times obtained with the Julia version of the IPOPT general solver for nonlinear programming [WB06] are also indicated. Both solvers were used with their respective default tolerances. With the same starting point, both solvers reach similar local solutions despite different stopping criteria.

We show the modeling of the problem with JuMP and its solving with IPOPT for n=1000.

cr_ipopt = Model(Ipopt.Optimizer); set_silent(cr_ipopt)
@variable(cr_ipopt, x[i=1:n], start = x0[i])

for k=1:n-2
@NLconstraint(cr_ipopt, 3x[k+1]^3 + 2x[k+2] - 5 + sin(x[k+1]-x[k+2])*sin(x[k+1]+x[k+2]) + 4x[k+1] - x[k]*exp(x[k]-x[k+1]) - 3 == 0)
end

@NLobjective(cr_ipopt, Min, sum(100(x[i]^2 - x[i+1])^2 + (x[i]-1)^2 for i=1:n-1))

optimize!(cr_ipopt)
termination_status(cr_ipopt)
LOCALLY_SOLVED::TerminationStatusCode = 4

The two solvers reach similar solutions.

Enlsip.sum_sq_residuals(cr_enlsip) ≈ objective_value(cr_ipopt)
true
Enlsip.solution(cr_enlsip) ≈ value.(cr_ipopt[:x])
true

The benchmark[BT] of the solving times (in seconds) for Enlsip.jl and Ipopt.jl is given in the following table:

Value of nEnlsip.jlIPOPT
103.616e-42.703e-3
1003.322e-24.759e-3
10002.325e02.520e-2
50003.172e21.490e-1
• BackendForwardDiff.jl https://juliadiff.org/ForwardDiff.jl/stable/
• HS80W. Hock and K. Schittkowski. Test Examples for Nonlinear Programming Codes, volume 187 of Lecture Notes in Economics and Mathematical Systems. Springer, second edition, 1980.
• LV99L. Lukšan and J. Vlček. Sparse and Partially Separable Test Problems for Unconstrained and Equality Constrained Optimization. Technical report 767, 1999.
• BTThe values shown in the table were obtained using the @btime macro from the BenchmarkTools.jl package.
• WB06A. Wächter and L. Biegler, On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming, Mathematical Programming, vol. 106(1), pages 25-57, 2006.