Performance

Preallocation

When calling our solver, Penelopt.jl starts with an allocation phase where all of the vectors, matrices, and sub-solvers needed to run the algorithm are built. The workspace is allocated in a struct called L2PenaltySolver. Our solver then proceeds to solve the provided nlp, in an allocation-free function named solve!.

In particular, this means that once a solver has been constructed, it can be reused across several calls to solve! without allocating any new memory, provided the problem dimensions ($n$ variables, $m$ constraints) stay the same. This is particularly useful:

when solving the same problem repeatedly from different starting points (e.g., in a multistart strategy);
when solving a sequence of closely related problems (e.g., inside a loop, for bilevel or parametric optimization);
when you want to avoid the garbage collector getting in the way of a tight benchmarking loop.

The convenience call

julia> stats = L2Penalty(nlp)

internally does exactly the following:

julia> solver = L2PenaltySolver(nlp)
julia> stats = PeneloptExecutionStats(nlp)
julia> solve!(solver, nlp, stats)

Separating these steps is the key to avoiding allocations: L2PenaltySolver(nlp) performs all the allocations up front (workspaces, factorization buffers, sub-solver states, etc.), and every subsequent call to solve! reuses that memory in place.

julia> using Penelopt

# Allocate the solver and the stats object once.
julia> solver = L2PenaltySolver(nlp)
julia> stats = PeneloptExecutionStats(nlp)

# Solve as many times as needed;
julia> solve!(solver, nlp, stats)
julia> solve!(solver, nlp, stats; x = x1)
julia> solve!(solver, nlp, stats; x = x2)

Resetting internal state

Some internal quantities of the solver (such as quasi-Newton approximations, the watchdog checkpoint, and factorization counters) persist between calls to solve!. Therefore, we strongly recommend calling SolverCore.reset! between calls to solve!.

julia> solve!(solver, nlp, stats)
julia> reset!(solver)
julia> solve!(solver, nlp, stats; x = x1)
julia> ...

Keyword arguments

Some options require allocations and are therefore passed to the solver structure construction call. The following options need to be passed to the solver constructor (passing them to solve! will cause failure).

r2n_m_monotone::Int = 12;
linear_solver::Sring = "ldlt".

For clarity, if you want to modify these options, while preallocating the workspace, you should do

julia> solver = L2PenaltySolver(nlp, r2n_m_monotone = 6, linear_solver = "mumps")
julia> ...

The `L2PenaltySolver` structure

The fields for the L2PenaltySolver are:

Field	Type	Description
`x`	`V`	current outer iterate $x_k$
`xn`	`V`	(internal) trial/next iterate buffer
`y`	`V`	current Lagrange multiplier estimate $y_k$
`cn`	`V`	(internal) buffer for the constraint values at a trial point
`dual_res`	`V`	buffer for the dual residual $\nabla f(x_k) + J(x_k)^Ty_k$
`s`	`V`	(internal) buffer for the computed step
`s0`	`V`	(internal) buffer used when computing least-squares multipliers
`∇fk`	`V`	buffer for $\nabla f(x_k)$
`temp_b`	`V`	(internal) scratch buffer of length $m$, used e.g. by the least-squares multiplier computation
`subsolver`	`PenaltyR2NSolver`	the R2N (inner-loop) solver structure used to (approximately) minimize the current penalized subproblem $f(x) + \tau_k\lVert c(x) \rVert_2$
`subpb`	`L2PenalizedProblem`	the penalized subproblem $f(x) + \tau_k\lVert c(x)\rVert_2$ associated with `nlp`
`substats`	`GenericExecutionStats`	the statistics object updated by `subsolver` at every R2N iteration

We mark some fields as internal which means that we discourage users to read/modify these in callbacks as they are only used for internal computations.

Here, T is the scalar type and V is the vector type of the problem (by default T == Float64 and V == Vector{Float64}; see the note on parametric types on the options page).

Nested solver structures

Because the algorithm is organized in nested loops (see the options terminology), solver.subsolver and solver.subsolver.subsolver expose the state of the inner loops.

solver.subsolver::PenaltyR2NSolver holds the state needed to minimize a single penalized subproblem $f(x) + \tau\lVert c(x)\rVert_2$ for a fixed $\tau$: the current inner iterate xk, its multiplier estimate y, the non-monotone objective history m_fh_hist, the watchdog checkpoint, and its own subsolver::MoreSorensenSolver and subpb::ShiftedL2PenalizedProblem.
solver.subsolver.subsolver::MoreSorensenSolver holds the buffers needed to compute a single step by (approximately) solving the trust-region subproblem via the Moré–Sorensen method: the right-hand side/solution buffers u1, u2, x1, x2, the assembled KKT-like matrix H, and the workspace used by the chosen linear solver (e.g., PenaltyLDLTWorkspace, PenaltyMUMPSWorkspace, ...).

All three levels can be constructed once (L2PenaltySolver, PenaltyR2NSolver, or MoreSorensenSolver, respectively) and reused across solves, following the same two-call pattern described above.

Fixed Variables

If some variables in your problem are fixed (i.e., some of your constraints are $x_i = c_i$ for some constants $c_i$), Penelopt.jl automatically removes these internally and solves a reduced problem.

using CUTEst, Penelopt

# NLP model with fixed variables
nlp = CUTEstModel("AIRCRFTA")

# You can check that a problem has fixed variables by accessing the `meta` of the NLP model
length(nlp.meta.ifix) > 0

true

The reformulation is done automatically when calling the solver, and the solution is then mapped back to the original problem. However, if you wish to use a preallocated solver you are responsible for constructing the reduced problem yourself.

The example below shows how to construct a preallocated solver for a problem with fixed variables.

# For the reformulation, it suffices to call the following function.
nlp_no_fixed = remove_fixed_variables(nlp)

# Now you can construct a preallocated solver for the reduced problem.
solver = L2PenaltySolver(nlp_no_fixed)
stats = PeneloptExecutionStats(nlp_no_fixed)

solve!(solver, nlp_no_fixed, stats)


# To retrieve the solution back you can do
solution_full = recover_full_solution(nlp_no_fixed, stats.solution)

8-element Vector{Float64}:
  0.005652720546060122
 -0.006537743733686216
 -0.0006212466687437924
 -0.1233355537445725
 -0.0003874221930345379
  0.1
  0.0
  0.0

Quasi-Newton Approximations

If the Hessian of the Lagrangian of your nonlinear programming problem is dense, ill-conditionned, expensive to compute, or inaccessible, you may be interested in replacing it with a quasi-Newton approximation.

Penelopt.jl offers you the possibility to run the optimization process with a Limited-memory BFGS approximation. You can simply pass a keyword argument when calling the solver (see the options page):

using CUTEst, Penelopt

# Construct your NLP model
nlp = CUTEstModel("HS6")

stats = L2Penalty(nlp; print_level = 1, qn_hessian_approximation = "bfgs")

┌ Info: 
│ This is Penelopt.jl v0.1.0.
│ Running with linear solver LDLFactorizations.jl v0.10.2.
│ 
│ CompactBFGSModel{Float64, Vector{Float64}, CUTEst.CUTEstModel{Float64}, NLPModels.NLPModelMeta{Float64, Vector{Float64}}, CompactBFGS{Float64, Vector{Float64}, Matrix{Float64}}} - A QuasiNewtonModel
│   Problem name: HS6
│    All variables: ████████████████████ 2      All constraints: ████████████████████ 1
│             free: ████████████████████ 2                 free: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
│            lower: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                lower: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
│            upper: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                upper: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
│          low/upp: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0              low/upp: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
│            fixed: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                fixed: ████████████████████ 1
│           infeas: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0               infeas: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
│             nnzh: ( 66.67% sparsity)   1               linear: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
│                                                     nonlinear: ████████████████████ 1
│                                                          nnzj: (  0.00% sparsity)   2
│                                                      lin_nnzj: (------% sparsity)
│                                                      nln_nnzj: (  0.00% sparsity)   2
│ 
│   Counters:
│              obj: ████████████████████ 1                 grad: ████████████████████ 1                 cons: ████████████████████ 1
│         cons_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0             cons_nln: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                 jcon: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
│            jgrad: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                  jac: ████████████████████ 1              jac_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
│          jac_nln: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                jprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0            jprod_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
│        jprod_nln: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0               jtprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0           jtprod_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
│       jtprod_nln: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                 hess: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                hprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
│            jhess: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0               jhprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
└ 
[ Info: ------------------------------------------------------------------------------------------------------
[ Info: Iter   sIter  Objective       pfeas       dfeas       τ           ptol        dtol        ‖x‖
[ Info: ------------------------------------------------------------------------------------------------------
[ Info: 0      0      +4.8400000e+00  8.15e+00    4.40e+00    1.00e+00    1.00e+00    8.15e-02    1.56e+00
[ Info: 1      9      +1.2912693e-03  7.35e-02    3.03e-02    1.00e+00    7.35e-04    3.03e-04    1.33e+00
[ Info: 2      2      +3.8532940e-08  1.85e-04    1.50e-04    1.00e+00    1.85e-06    1.50e-06    1.41e+00
[ Info: 3      2      +1.8634535e-16  8.32e-11    5.06e-08    1.00e+00    8.05e-08    1.36e-07    1.41e+00
┌ Info: 
│ Number of Iterations: 3
│ 
│ 
│ Objective...........: +1.863453482150732e-16
│ Primal Feasibility..:  8.323786104824649e-11
│ Dual Feasibility....:  5.055633020672882e-08
│ 
│ 
└ EXIT: first_order.

Default options

Some default options have a different value when using our solver with a quasi-Newton approximation. Those options are

r2n_η2::T = 0.1 which default to 0.9 with a quasi-Newton approximation.

L-BFGS options

You can pass the following additional keyword arguments to customize your approximation. Those are

qn_mem::Int = 6: memory parameter for the limited-memory approximation.
qn_scaling::Bool = true: whether we scale $B_0 = \gamma I$ with $\gamma = y^Ty / s^T y$, where $s$ is the step and $y$ is the difference of the two last gradients of the Lagrangian.
qn_max_skip::Int = 2 (advanced): if we skipped a pair $(s, y)$ more than max_skip times in a row, we reset the approximation.

If you wish to use a quasi-Newton approximation together with the preallocation feature you should do the following instead:

julia> using CUTEst, Penelopt

# Construct your NLP model
julia> nlp = CUTEstModel("HS6")

# Construct a "quasi-Newton" NLP structure first
julia> nlp_bfgs = CompactBFGSModel(nlp; mem = 6, scaling = true, max_skip = 2)

# The solver will automatically use the quasi-Newton approximation now
julia> solver = L2PenaltySolver(nlp_bfgs)
julia> stats = PeneloptExecutionStats(nlp_bfgs)

# Solve as many times as needed;
julia> solve!(solver, nlp, stats)

Quasi-Newton approximations and fixed variables

If your problem has fixed variables, you should first remove them before constructing a quasi-Newton approximation. This is because the quasi-Newton approximation is built for the reduced problem and it will not be valid for the original problem.