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].

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).

solve!(model{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:

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

solution(model)

status(model)

constraints_values(model)

sum_sq_residuals(model)

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.321
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: