Optimization

Paul Schrimpf

2026-03-09

Optimization

Optimization

\[ \max_{x} f(x) \]

Overview of Algorithms

• https://schrimpf.github.io/AnimatedOptimization.jl/optimization/

Example Problem

using Statistics, LinearAlgebra

function cueiv(β; n=1000, σ=0.1, γ=[I; ones(2,length(β))], ρ=0.5)
  z = randn(n, size(γ)[1])
  endo = randn(n, length(β))
  x = z*γ .+ endo
  ϵ = σ*(randn(n)*sqrt(1.0-ρ^2).+endo[:,1]*ρ)
  y =  x*β .+ ϵ
  g(β) = (y - x*β).*z
  function cueobj(β)
    G = g(β)
    Eg=mean(G,dims=1)
    W = inv(cov(G))
    (n*Eg*W*Eg')[1]
  end
  return(cueobj, x, y, z)
end

cueiv (generic function with 1 method)

using Plots
β = ones(2)
f, x, y ,z = cueiv(β;σ=0.1,ρ=0.9)
b1 = range(0,2,length=100)
b2 = range(0,2,length=100)
fval = ((x,y)->f([x,y])).(b1,b2')
Plots.plotly()
contourf(b1,b2,fval)

Packages

We focus on algorithms useful for small to medium scale smooth nonlinear problems. Variants of Newton’s method will generally be the best choice for such problems.

Optim.jl

Pure Julia implementation of standard set of optimization algorithms. This is generally a good starting point, especially for unconstrained nonlinear problems without too many variables (maybe up to a few hundred).

• Optim docs

import Optim, ADTypes, ForwardDiff
using PrettyTables
β = ones(2)
f = cueiv(β;σ=0.1,ρ=0.5)[1]
β0 = zero(β)

sol=Optim.optimize(f, β0, Optim.NelderMead())
sol=Optim.optimize(f,β0, Optim.LBFGS(), autodiff=ADTypes.AutoForwardDiff())
sol=Optim.optimize(f,β0, Optim.LBFGS(m=20, linesearch=Optim.LineSearches.BackTracking()), autodiff=ADTypes.AutoForwardDiff())
sol=Optim.optimize(f,β0,Optim.NewtonTrustRegion(), autodiff=ADTypes.AutoForwardDiff())

 * Status: success

 * Candidate solution
    Final objective value:     6.484982e-02

 * Found with
    Algorithm:     Newton's Method (Trust Region)

 * Convergence measures
    |x - x'|               = 1.44e-09 ≰ 0.0e+00
    |x - x'|/|x'|          = 1.44e-09 ≰ 0.0e+00
    |f(x) - f(x')|         = 3.93e-13 ≰ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 6.06e-12 ≰ 0.0e+00
    |g(x)|                 = 1.84e-11 ≤ 1.0e-08

 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    21
    f(x) calls:    22
    ∇f(x) calls:   22
    ∇f(x)ᵀv calls: 0
    ∇²f(x) calls:  15
    ∇²f(x)v calls: 0

Tougher in Higher Dimension

β = ones(20)
f =  cueiv(β)[1]
β0 = zero(β)
sol=Optim.optimize(f, β0, Optim.NelderMead())
sol=Optim.optimize(f,β0, Optim.LBFGS(), autodiff=ADTypes.AutoForwardDiff())
sol=Optim.optimize(f,β0, Optim.LBFGS(m=20, linesearch=Optim.LineSearches.BackTracking()), autodiff=ADTypes.AutoForwardDiff())
sol=Optim.optimize(f,β0,Optim.NewtonTrustRegion(), autodiff=ADTypes.AutoForwardDiff())

 * Status: success

 * Candidate solution
    Final objective value:     1.204639e+00

 * Found with
    Algorithm:     Newton's Method (Trust Region)

 * Convergence measures
    |x - x'|               = 2.18e-11 ≰ 0.0e+00
    |x - x'|/|x'|          = 2.17e-11 ≰ 0.0e+00
    |f(x) - f(x')|         = 1.62e-14 ≰ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 1.35e-14 ≰ 0.0e+00
    |g(x)|                 = 1.66e-10 ≤ 1.0e-08

 * Work counters
    Seconds run:   3  (vs limit Inf)
    Iterations:    33
    f(x) calls:    34
    ∇f(x) calls:   34
    ∇f(x)ᵀv calls: 0
    ∇²f(x) calls:  25
    ∇²f(x)v calls: 0

Optimization.jl

A package that provides a uniform interface to many other optimization packages. A good choice for easy experimentation and trying different packages.

• Optimization docs
• Unified interface to many packages (like NonlinearSolve.jl, but less polished)

JuliaSmoothOptimizers

A smaller set of developers and users than Optim, but has some good ideas and algorithms. Multiple options for nonlinear constrained optimization, some with promising benchmarks.. Design makes using specialized linear solvers inside nonlinear algorithms convenient.

• Organization README

using ADNLPModels
import DCISolver, JSOSolvers, Percival
β = ones(2)
f =  cueiv(β; γ = I + zeros(2,2))[1]
β0 = zero(β)

nlp = ADNLPModel(f, β0)
stats=JSOSolvers.lbfgs(nlp)
print(stats)

stats=JSOSolvers.tron(nlp)
print(stats)

ecnlp = ADNLPModel(f,β0, b->(length(β0) - sum(b)), zeros(1), zeros(1))
stats=DCISolver.dci(ecnlp)
print(stats)

icnlp = ADNLPModel(f,β0, b->(length(β0) - sum(b)), -ones(1), ones(1))
stats=Percival.percival(icnlp)
print(stats)

JuMP

• JuMP is a modelling language for optimization in Julia
• many solvers
• best with problem written in special JuMP syntax
  • efficiently calculate derivatives and recognize special problem structure such as linearity, sparsity, etc.

using JuMP
import MadNLP, Ipopt
k = 10
β = ones(k)
f,x,y,z =  cueiv(β; γ = I(k))
β0 = zero(β)
n,k = size(x)

# sub-optimal usage of JuMP
m = Model(()->MadNLP.Optimizer(print_level=MadNLP.INFO, max_iter=100))
#m = Model(Ipopt.Optimizer)
@variable(m, β[1:k])
@operator(m, op_f, k, (x...)->f(vcat(x...))) # hides internals of f from JuMP
@objective(m, Min, op_f(β...))
optimize!(m)
JuMP.value.(β)

This is MadNLP version v0.9.1, running with MUMPS v5.8.2

Number of nonzeros in constraint Jacobian............:        0
Number of nonzeros in Lagrangian Hessian.............:        0

Total number of variables............................:       10
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:        0
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du inf_compl lg(mu) lg(rg) alpha_pr ir ls
   0  2.9567132e+02 0.00e+00 5.08e+01 0.00e+00  -1.0     -   0.00e+00  1  0 
   1  2.9336479e+02 0.00e+00 1.34e+01 0.00e+00  -1.0     -   7.81e-03  1  8f
   2  2.7430036e+02 0.00e+00 7.60e+00 0.00e+00  -1.0     -   1.00e+00  1  1f
   3  2.6900612e+02 0.00e+00 1.41e+01 0.00e+00  -1.0     -   1.00e+00  1  1f
   4  2.4690168e+02 0.00e+00 4.75e+00 0.00e+00  -1.0     -   1.00e+00  1  1f
   5  2.3974535e+02 0.00e+00 3.70e+00 0.00e+00  -1.0     -   1.00e+00  1  1f
   6  2.3812993e+02 0.00e+00 1.72e+00 0.00e+00  -1.0     -   1.00e+00  1  1f
   7  2.3742262e+02 0.00e+00 3.62e-01 0.00e+00  -1.0     -   1.00e+00  1  1f
   8  2.3739492e+02 0.00e+00 9.98e-02 0.00e+00  -1.7     -   1.00e+00  1  1h
   9  2.3739198e+02 0.00e+00 2.00e-02 0.00e+00  -2.5     -   1.00e+00  1  1h
iter    objective    inf_pr   inf_du inf_compl lg(mu) lg(rg) alpha_pr ir ls
  10  2.3739189e+02 0.00e+00 7.27e-03 0.00e+00  -3.8     -   1.00e+00  1  1h
  11  2.3739188e+02 0.00e+00 4.61e-03 0.00e+00  -3.8     -   1.00e+00  1  1h
  12  2.3739185e+02 0.00e+00 3.10e-03 0.00e+00  -3.8     -   1.00e+00  1  1h
  13  2.3739182e+02 0.00e+00 1.52e-02 0.00e+00  -3.8     -   1.25e-01  1  4h
  14  2.3739177e+02 0.00e+00 2.18e-02 0.00e+00  -3.8     -   1.25e-01  1  4h
  15  2.3739165e+02 0.00e+00 5.16e-02 0.00e+00  -3.8     -   2.50e-01  1  3h
  16  2.3739162e+02 0.00e+00 6.56e-02 0.00e+00  -3.8     -   6.25e-02  1  5h
  17  2.3739105e+02 0.00e+00 6.81e-02 0.00e+00  -3.8     -   3.12e-02  1  6h
  18  2.3739084e+02 0.00e+00 9.20e-02 0.00e+00  -3.8     -   1.56e-02  1  7h
  19  2.3737958e+02 0.00e+00 2.53e-01 0.00e+00  -3.8     -   3.12e-02  1  6h
iter    objective    inf_pr   inf_du inf_compl lg(mu) lg(rg) alpha_pr ir ls
  20  2.3737398e+02 0.00e+00 6.52e-01 0.00e+00  -3.8     -   7.63e-06  1 18f
  21  2.3737397e+02 0.00e+00 6.65e-01 0.00e+00  -3.8     -   2.44e-04  1 13h
  22  2.3737141e+02 0.00e+00 1.14e+00 0.00e+00  -3.8     -   6.25e-02  1  5h
  23  2.3736885e+02 0.00e+00 2.08e+00 0.00e+00  -3.8     -   9.77e-04  1 11f
  24  2.3736828e+02 0.00e+00 2.90e+00 0.00e+00  -3.8     -   4.88e-04  1 12f
  25  2.3736820e+02 0.00e+00 3.00e+00 0.00e+00  -3.8     -   4.88e-04  1 12h
  26  2.3736794e+02 0.00e+00 3.09e+00 0.00e+00  -3.8     -   3.91e-03  1  9h
  27  2.0262079e+02 0.00e+00 1.36e+03 0.00e+00  -3.8     -   9.77e-04  1 11f
  28  2.0038824e+02 0.00e+00 1.58e+03 0.00e+00  -3.8     -   7.81e-03  1  8f
  29  1.9276671e+02 0.00e+00 1.94e+03 0.00e+00  -3.8     -   1.95e-03  1 10f
iter    objective    inf_pr   inf_du inf_compl lg(mu) lg(rg) alpha_pr ir ls
  30  1.8786417e+02 0.00e+00 2.02e+03 0.00e+00  -3.8     -   3.91e-03  1  9f
  31  1.8725759e+02 0.00e+00 1.98e+03 0.00e+00  -3.8     -   7.81e-03  1  8f
  32  1.0287665e+02 0.00e+00 2.08e+03 0.00e+00  -3.8     -   1.00e+00  1  1f
  33  5.7571027e+01 0.00e+00 1.77e+03 0.00e+00  -3.8     -   1.00e+00  1  1f
  34  5.4227201e+01 0.00e+00 1.76e+03 0.00e+00  -3.8     -   1.00e+00  1  1f
  35  1.8100694e+01 0.00e+00 1.28e+03 0.00e+00  -3.8     -   5.00e-01  1  2f
  36  1.7081915e+01 0.00e+00 1.47e+03 0.00e+00  -3.8     -   1.00e+00  1  1f
  37  8.4050940e-01 0.00e+00 4.25e+02 0.00e+00  -3.8     -   1.00e+00  1  1f
  38  4.5789992e-01 0.00e+00 3.04e+02 0.00e+00  -3.8     -   1.00e+00  1  1f
  39  8.4414662e-03 0.00e+00 3.51e+01 0.00e+00  -3.8     -   1.00e+00  1  1h
iter    objective    inf_pr   inf_du inf_compl lg(mu) lg(rg) alpha_pr ir ls
  40  6.8214059e-04 0.00e+00 1.07e+01 0.00e+00  -3.8     -   1.00e+00  1  1h
  41  3.4925261e-05 0.00e+00 2.37e+00 0.00e+00  -3.8     -   1.00e+00  1  1h
  42  1.0977887e-06 0.00e+00 4.23e-01 0.00e+00  -3.8     -   1.00e+00  1  1h
  43  1.2289961e-08 0.00e+00 4.16e-02 0.00e+00  -3.8     -   1.00e+00  1  1h
  44  7.6589841e-11 0.00e+00 3.03e-03 0.00e+00  -3.8     -   1.00e+00  1  1h
  45  7.0288306e-13 0.00e+00 2.46e-04 0.00e+00  -3.8     -   1.00e+00  1  1h
  46  6.6725883e-15 0.00e+00 3.41e-05 0.00e+00  -5.7     -   1.00e+00  1  1h
  47  4.5966786e-16 0.00e+00 9.75e-06 0.00e+00  -5.7     -   1.00e+00  1  1h
  48  1.0202229e-18 0.00e+00 3.48e-07 0.00e+00  -8.6     -   1.00e+00  1  1h
  49  7.5939896e-20 0.00e+00 7.80e-08 0.00e+00  -8.6     -   1.00e+00  1  1h
iter    objective    inf_pr   inf_du inf_compl lg(mu) lg(rg) alpha_pr ir ls
  50  1.9915077e-21 0.00e+00 1.64e-08 0.00e+00  -8.6     -   1.00e+00  1  1h
  51  1.9075041e-22 0.00e+00 5.72e-09 0.00e+00  -9.0     -   1.00e+00  1  1h

Number of Iterations....: 51

                                   (scaled)                 (unscaled)
Objective...............:   1.9075041404754786e-22    1.9075041404754786e-22
Dual infeasibility......:   5.7216939700066386e-09    5.7216939700066386e-09
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   5.7216939700066386e-09    5.7216939700066386e-09

Number of objective function evaluations              = 286
Number of objective gradient evaluations              = 52
Number of constraint evaluations                      = 286
Number of constraint Jacobian evaluations             = 1
Number of Lagrangian Hessian evaluations              = 0
Number of KKT factorizations                          = 52
Number of KKT backsolves                              = 136

Total wall secs in initialization                     =  0.635 s
Total wall secs in linear solver                      =  0.012 s
Total wall secs in NLP function evaluations           =  0.172 s
Total wall secs in solver (w/o init./fun./lin. alg.)  =  5.136 s
Total wall secs                                       =  5.955 s

EXIT: Optimal Solution Found (tol = 1.0e-08).

10-element Vector{Float64}:
 1.0003842241777336
 1.0000774341695478
 1.000784391027142
 1.0034346882976357
 1.0000664938514874
 0.9937734717435119
 0.9996355470883188
 1.0027586733191065
 0.9999785590469339
 0.9999614128181697

# better usage of JuMP
m2 = Model(()->MadNLP.Optimizer(print_level=MadNLP.INFO, max_iter=100))
#m2 = Model(Ipopt.Optimizer)
#set_attribute(m2, "print_level", 5)
n, k = size(x)
@variable(m2, β2[1:k])
g = (y - x*β2).*z
Eg=mean(g,dims=1)
invW = cov(g)
@variable(m2, Wg[1:k])
@constraint(m2, invW*Wg .== Eg')
@objective(m2, Min, n*dot(Eg,Wg))
optimize!(m2)
JuMP.value.(β2)

This is MadNLP version v0.9.1, running with MUMPS v5.8.2

Number of nonzeros in constraint Jacobian............:      200
Number of nonzeros in Lagrangian Hessian.............:     1750

Total number of variables............................:       20
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:       10
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du inf_compl lg(mu) lg(rg) alpha_pr ir ls
   0  0.0000000e+00 1.19e+00 4.01e+00 0.00e+00  -1.0     -   0.00e+00  1  0 
   1  2.9156891e+02 1.26e-02 4.39e+00 0.00e+00  -1.0    2.0  1.00e+00  1  1h
   2  2.8127669e+02 3.30e-03 3.73e+00 0.00e+00  -1.0    1.5  1.00e+00  1  1h
   3  2.6384809e+02 1.84e-02 2.83e+00 0.00e+00  -1.0    1.0  1.00e+00  1  1h
   4  2.4935484e+02 6.37e-02 5.32e+00 0.00e+00  -1.0    0.6  1.00e+00  1  1h
   5  2.4162040e+02 4.94e-02 4.14e+00 0.00e+00  -1.0    0.1  1.00e+00  1  1h
   6  2.3631989e+02 1.37e-02 1.15e+00 0.00e+00  -1.0     -   1.00e+00  1  1h
   7  2.3716692e+02 1.50e-04 1.84e-02 0.00e+00  -1.0   -0.4  1.00e+00  1  1h
   8  2.3682077e+02 2.18e-03 1.82e-01 0.00e+00  -3.8   -0.9  1.00e+00  1  1h
   9  1.6554877e+02 2.07e-01 1.73e+01 0.00e+00  -3.8   -1.3  1.00e+00  1  1h
iter    objective    inf_pr   inf_du inf_compl lg(mu) lg(rg) alpha_pr ir ls
  10  2.2830477e+02 2.24e-02 1.01e+00 0.00e+00  -3.8    1.8  1.00e+00  1  1h
  11  2.3639656e+02 1.46e-04 1.38e-01 0.00e+00  -3.8    1.3  1.00e+00  1  1h
  12  2.3632971e+02 2.32e-04 1.03e-01 0.00e+00  -3.8    0.8  1.00e+00  1  1h
  13  2.3508250e+02 2.98e-03 2.66e-01 0.00e+00  -3.8    0.4  1.00e+00  1  1h
  14  2.3586954e+02 8.18e-04 1.68e-01 0.00e+00  -3.8    0.8  1.00e+00  1  1h
  15  2.3099041e+02 9.96e-03 9.10e-01 0.00e+00  -3.8    0.3  1.00e+00  1  1h
  16  2.2887104e+02 1.12e-02 1.36e+00 0.00e+00  -3.8    0.7  1.00e+00  1  1h
  17  2.3378155e+02 2.15e-03 4.57e-01 0.00e+00  -3.8    1.2  1.00e+00  1  1h
  18  2.2155945e+02 1.54e-02 3.34e+00 0.00e+00  -3.8    0.7  1.00e+00  1  1h
  19  1.4161942e+02 3.96e-02 3.47e+01 0.00e+00  -3.8    1.1  1.00e+00  1  1h
iter    objective    inf_pr   inf_du inf_compl lg(mu) lg(rg) alpha_pr ir ls
  20 -2.5025059e+03 1.25e+01 1.16e+03 0.00e+00  -3.8    2.4  1.00e+00  1  1H
  21 -1.0947655e+03 3.75e+00 1.99e+02 0.00e+00  -3.8    2.9  1.00e+00  1  1h
  22 -3.6353218e+02 1.21e+00 6.79e+01 0.00e+00  -3.8     -   1.00e+00  1  1h
  23 -3.3038431e+02 6.33e-01 5.28e+01 0.00e+00  -3.8     -   1.00e+00  1  1h
  24 -1.8294657e+02 2.24e-01 4.03e+01 0.00e+00  -3.8    2.4  1.00e+00  1  1h
  25 -1.0672755e+02 7.98e-02 3.65e+01 0.00e+00  -3.8     -   5.00e-01  1  2h
  26 -1.8536937e+01 1.82e-02 1.67e+01 0.00e+00  -3.8     -   1.00e+00  1  1h
  27  1.5653516e-01 8.04e-04 1.65e+01 0.00e+00  -3.8     -   1.00e+00  1  1h
  28  7.6419111e-05 1.02e-05 1.70e-01 0.00e+00  -3.8     -   1.00e+00  1  1h
  29 -1.3393526e-12 5.00e-09 1.33e-05 0.00e+00  -3.8     -   1.00e+00  1  1h
iter    objective    inf_pr   inf_du inf_compl lg(mu) lg(rg) alpha_pr ir ls
  30  0.0000000e+00 3.10e-16 3.92e-12 0.00e+00  -8.6     -   1.00e+00  1  1h

Number of Iterations....: 30

                                   (scaled)                 (unscaled)
Objective...............:   0.0000000000000000e+00    0.0000000000000000e+00
Dual infeasibility......:   3.9171273809825470e-12    4.6662054384502229e-11
Constraint violation....:   3.1034825007030597e-16    3.1034825007030597e-16
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   3.2883007662730172e-13    3.9171273809825470e-12

Number of objective function evaluations              = 33
Number of objective gradient evaluations              = 31
Number of constraint evaluations                      = 33
Number of constraint Jacobian evaluations             = 31
Number of Lagrangian Hessian evaluations              = 30
Number of KKT factorizations                          = 66
Number of KKT backsolves                              = 32

Total wall secs in initialization                     =  0.010 s
Total wall secs in linear solver                      =  0.001 s
Total wall secs in NLP function evaluations           =  1.600 s
Total wall secs in solver (w/o init./fun./lin. alg.)  =  1.614 s
Total wall secs                                       =  3.224 s

EXIT: Optimal Solution Found (tol = 1.0e-08).

10-element Vector{Float64}:
 1.0003842241777279
 1.0000774341695595
 1.0007843910271286
 1.0034346882976517
 1.0000664938514778
 0.9937734717435205
 0.9996355470883034
 1.002758673319113
 0.9999785590469424
 0.9999614128181896

# or Empirical Likelihood
m3 = Model(()->MadNLP.Optimizer(print_level=MadNLP.INFO, max_iter=100))
@variable(m3, β3[1:k])
@variable(m3, 1e-8 <= p[1:n] <= 1-1e-8)
@constraint(m3, sum(p) == 1)
@objective(m3, Max, 1/n*sum(log.(p)))
JuMP.set_start_value.(p, 1/n)
g3 = (y - x*β3).*z
@constraint(m3, p'*g3 .== 0)
optimize!(m3)
JuMP.value.(β3)

This is MadNLP version v0.9.1, running with MUMPS v5.8.2

Number of nonzeros in constraint Jacobian............:   211000
Number of nonzeros in Lagrangian Hessian.............:   101000

Total number of variables............................:     1010
                     variables with only lower bounds:        0
                variables with lower and upper bounds:     1000
                     variables with only upper bounds:        0
Total number of equality constraints.................:       11
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du inf_compl lg(mu) lg(rg) alpha_pr ir ls
   0  4.6051702e+00 1.19e+01 1.39e-16 9.90e-01  -1.0     -   0.00e+00  1  0 
   1  6.9077553e+00 1.07e+00 8.10e-01 9.18e-02  -1.0     -   1.00e+00  1  1h
   2  6.9077553e+00 2.81e-13 7.11e-14 1.00e-01  -1.0     -   1.00e+00  1  1h
   3  6.9077553e+00 2.51e-14 4.84e-14 2.83e-03  -2.5     -   1.00e+00  1  1h
   4  6.9077553e+00 1.40e-14 2.66e-15 1.84e-06  -5.7     -   1.00e+00  1  1h
   5  6.9077553e+00 7.55e-15 4.44e-16 9.09e-10  -9.0     -   1.00e+00  1  1h

Number of Iterations....: 5

                                   (scaled)                 (unscaled)
Objective...............:   6.9077552789822052e+00    6.9077552789822052e+00
Dual infeasibility......:   4.4408920985006262e-16    4.4408920985006262e-16
Constraint violation....:   7.5495165674510645e-15    7.5495165674510645e-15
Complementarity.........:   9.0909090949664109e-10    9.0909090949664109e-10
Overall NLP error.......:   9.0909090949664109e-10    9.0909090949664109e-10

Number of objective function evaluations              = 6
Number of objective gradient evaluations              = 6
Number of constraint evaluations                      = 6
Number of constraint Jacobian evaluations             = 6
Number of Lagrangian Hessian evaluations              = 5
Number of KKT factorizations                          = 6
Number of KKT backsolves                              = 6

Total wall secs in initialization                     =  0.018 s
Total wall secs in linear solver                      =  0.007 s
Total wall secs in NLP function evaluations           =  0.002 s
Total wall secs in solver (w/o init./fun./lin. alg.)  =  0.005 s
Total wall secs                                       =  0.031 s

EXIT: Optimal Solution Found (tol = 1.0e-08).

10-element Vector{Float64}:
 1.0003842241777274
 1.0000774341695535
 1.00078439102713
 1.003434688297651
 1.0000664938514794
 0.9937734717435208
 0.9996355470883086
 1.002758673319118
 0.9999785590469416
 0.9999614128181843

ExaModels

• ExaModels.jl
• Somewhat similar to JuMP, but focused on nonlinear programs
  • less friendly syntax
  • better preformance for some problems
• GPU support

Ipopt

An open source solver written that works well on large constrained nonlinear problems. Can be used directly or through JuMP, Optimization.jl or JuliaSmoothOptimizers interfaces.

MadNLP

Interior point solver written in Julia.

Knitro

A commercial solver with state of the art performance on many problems. Maybe worth the cost for large constrained nonlinear problems. Can be used directly or through JuMP, Optimization.jl or JuliaSmoothOptimizers interfaces.

Optimisers

Optimisers.jl has a collection of gradient descent variations designed for use in machine learning. A reasonable choice for problems with a very large number of variables.

Covariance Matrix Adaption Evolution Strategy

CMAES can be effective for problems with many local minimum while retaining some of the benefits of Newton’s method. It can even deal with noisy and discontinuous objective functions. CMAEvolutionStrategy is a Julia implementation. GCMAES implements both the original CMAES algorithm and a variation that also uses gradient information.

import CMAEvolutionStrategy, GCMAES
β = ones(10)
f =  cueiv(β; γ = I + zeros(length(β),length(β)))[1]
β0 = zero(β)

out = CMAEvolutionStrategy.minimize(f, β0, 1.0;
                                    lower = nothing,
                                    upper = nothing,
                                    popsize = 4 + floor(Int, 3*log(length(β0))),
                                    verbosity = 1,
                                    seed = rand(UInt),
                                    maxtime = nothing,
                                    maxiter = 500,
                                    maxfevals = nothing,
                                    ftarget = nothing,
                                    xtol = nothing,
                                    ftol = 1e-11)
@show CMAEvolutionStrategy.xbest(out)

xmin, fmin, status = GCMAES.minimize(f, β0, 1.0, fill(-10.,length(β0)), fill(10., length(β0)), maxiter=500)
@show xmin

import ForwardDiff
∇f(β) = ForwardDiff.gradient(f, β)
xmin, fmin, status=GCMAES.minimize((f, ∇f), β0, 1.0, fill(-10.,length(β0)), fill(10., length(β0)), maxiter=500)
@show xmin