The nonlinear constrained optimization interface in
Optim
assumes that the user can write the optimization
problem in the following way.
For equality constraints on $x_j$ or $c(x)_j$ you set those particular entries of bounds to be equal, $l_j=u_j$. Likewise, setting $l_j=-\infty$ or $u_j=\infty$ means that the constraint is unbounded from below or above respectively.
using Optim, NLSolversBase #hide
import NLSolversBase: clear! #hide
IPNewton
¶We will go through examples on how to use the constraints interface with the interior-point Newton optimization algorithm IPNewton.
Throughout these examples we work with the standard Rosenbrock function. The objective and its derivatives are given by
fun(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
function fun_grad!(g, x)
g[1] = -2.0 * (1.0 - x[1]) - 400.0 * (x[2] - x[1]^2) * x[1]
g[2] = 200.0 * (x[2] - x[1]^2)
end
function fun_hess!(h, x)
h[1, 1] = 2.0 - 400.0 * x[2] + 1200.0 * x[1]^2
h[1, 2] = -400.0 * x[1]
h[2, 1] = -400.0 * x[1]
h[2, 2] = 200.0
end;
To solve a constrained optimization problem we call the optimize
method
optimize(d::AbstractObjective, constraints::AbstractConstraints, initial_x::Tx, method::ConstrainedOptimizer, options::Options)
We can create instances of AbstractObjective
and
AbstractConstraints
using the types TwiceDifferentiable
and
TwiceDifferentiableConstraints
from the package NLSolversBase.jl
.
We want to optimize the Rosenbrock function in the box
$-0.5 \leq x \leq 0.5$, starting from the point $x_0=(0,0)$.
Box constraints are defined using, for example,
TwiceDifferentiableConstraints(lx, ux)
.
x0 = [0.0, 0.0]
df = TwiceDifferentiable(fun, fun_grad!, fun_hess!, x0)
lx = [-0.5, -0.5]; ux = [0.5, 0.5]
dfc = TwiceDifferentiableConstraints(lx, ux)
res = optimize(df, dfc, x0, IPNewton())
* Status: success * Candidate solution Final objective value: 2.500000e-01 * Found with Algorithm: Interior Point Newton * Convergence measures |x - x'| = 4.39e-10 ≰ 0.0e+00 |x - x'|/|x'| = 8.79e-10 ≰ 0.0e+00 |f(x) - f(x')| = 0.00e+00 ≤ 0.0e+00 |f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00 |g(x)| = 1.00e+00 ≰ 1.0e-08 * Work counters Seconds run: 1 (vs limit Inf) Iterations: 43 f(x) calls: 68 ∇f(x) calls: 68
Like the rest of Optim, you can also use autodiff=:forward
and just pass in
fun
.
If we only want to set lower bounds, use ux = fill(Inf, 2)
ux = fill(Inf, 2)
dfc = TwiceDifferentiableConstraints(lx, ux)
clear!(df)
res = optimize(df, dfc, x0, IPNewton())
* Status: success (objective increased between iterations) * Candidate solution Final objective value: 7.987239e-20 * Found with Algorithm: Interior Point Newton * Convergence measures |x - x'| = 3.54e-10 ≰ 0.0e+00 |x - x'|/|x'| = 3.54e-10 ≰ 0.0e+00 |f(x) - f(x')| = 2.40e-19 ≰ 0.0e+00 |f(x) - f(x')|/|f(x')| = 3.00e+00 ≰ 0.0e+00 |g(x)| = 8.83e-09 ≤ 1.0e-08 * Work counters Seconds run: 0 (vs limit Inf) Iterations: 35 f(x) calls: 63 ∇f(x) calls: 63
An unconstrained problem can be defined either by passing
Inf
bounds or empty arrays.
Note that we must pass the correct type information to the empty lx
and ux
lx = fill(-Inf, 2); ux = fill(Inf, 2)
dfc = TwiceDifferentiableConstraints(lx, ux)
clear!(df)
res = optimize(df, dfc, x0, IPNewton())
lx = Float64[]; ux = Float64[]
dfc = TwiceDifferentiableConstraints(lx, ux)
clear!(df)
res = optimize(df, dfc, x0, IPNewton())
* Status: success * Candidate solution Final objective value: 5.998937e-19 * Found with Algorithm: Interior Point Newton * Convergence measures |x - x'| = 1.50e-09 ≰ 0.0e+00 |x - x'|/|x'| = 1.50e-09 ≰ 0.0e+00 |f(x) - f(x')| = 1.80e-18 ≰ 0.0e+00 |f(x) - f(x')|/|f(x')| = 3.00e+00 ≰ 0.0e+00 |g(x)| = 7.92e-09 ≤ 1.0e-08 * Work counters Seconds run: 0 (vs limit Inf) Iterations: 34 f(x) calls: 63 ∇f(x) calls: 63
We now consider the Rosenbrock problem with a constraint on $$ c(x)_1 = x_1^2 + x_2^2. $$
We pass the information about the constraints to optimize
by defining a vector function c(x)
and its Jacobian J(x)
.
The Hessian information is treated differently, by considering the
Lagrangian of the corresponding slack-variable transformed
optimization problem. This is similar to how the CUTEst
library works.
Let $H_j(x)$ represent the Hessian of the $j$th component
$c(x)_j$ of the generic constraints.
and $\lambda_j$ the corresponding dual variable in the
Lagrangian. Then we want the constraint
object to
add the values of $H_j(x)$ to the Hessian of the objective,
weighted by $\lambda_j$.
The Julian form for the supplied function $c(x)$ and the derivative information is then added in the following way.
con_c!(c, x) = (c[1] = x[1]^2 + x[2]^2; c)
function con_jacobian!(J, x)
J[1,1] = 2*x[1]
J[1,2] = 2*x[2]
J
end
function con_h!(h, x, λ)
h[1,1] += λ[1]*2
h[2,2] += λ[1]*2
end;
Note that con_h!
adds the λ
-weighted Hessian value of each
element of c(x)
to the Hessian of fun
.
We can then optimize the Rosenbrock function inside the ball of radius $0.5$.
lx = Float64[]; ux = Float64[]
lc = [-Inf]; uc = [0.5^2]
dfc = TwiceDifferentiableConstraints(con_c!, con_jacobian!, con_h!,
lx, ux, lc, uc)
res = optimize(df, dfc, x0, IPNewton())
* Status: success * Candidate solution Final objective value: 2.966216e-01 * Found with Algorithm: Interior Point Newton * Convergence measures |x - x'| = 0.00e+00 ≤ 0.0e+00 |x - x'|/|x'| = 0.00e+00 ≤ 0.0e+00 |f(x) - f(x')| = 0.00e+00 ≤ 0.0e+00 |f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00 |g(x)| = 7.71e-01 ≰ 1.0e-08 * Work counters Seconds run: 0 (vs limit Inf) Iterations: 28 f(x) calls: 109 ∇f(x) calls: 109
We can add a lower bound on the constraint, and thus optimize the objective on the annulus with inner and outer radii $0.1$ and $0.5$ respectively.
lc = [0.1^2]
dfc = TwiceDifferentiableConstraints(con_c!, con_jacobian!, con_h!,
lx, ux, lc, uc)
res = optimize(df, dfc, x0, IPNewton())
┌ Warning: Initial guess is not an interior point └ @ Optim ~/work/Optim.jl/Optim.jl/src/multivariate/solvers/constrained/ipnewton/ipnewton.jl:116 Stacktrace: [1] initial_state(method::IPNewton{typeof(Optim.backtrack_constrained_grad), Symbol}, options::Optim.Options{Float64, Nothing}, d::TwiceDifferentiable{Float64, Vector{Float64}, Matrix{Float64}, Vector{Float64}}, constraints::TwiceDifferentiableConstraints{typeof(Main.var"##226".con_c!), typeof(Main.var"##226".con_jacobian!), typeof(Main.var"##226".con_h!), Float64}, initial_x::Vector{Float64}) @ Optim ~/work/Optim.jl/Optim.jl/src/multivariate/solvers/constrained/ipnewton/ipnewton.jl:117 [2] optimize @ ~/work/Optim.jl/Optim.jl/src/multivariate/solvers/constrained/ipnewton/interior.jl:229 [inlined] [3] optimize(d::TwiceDifferentiable{Float64, Vector{Float64}, Matrix{Float64}, Vector{Float64}}, constraints::TwiceDifferentiableConstraints{typeof(Main.var"##226".con_c!), typeof(Main.var"##226".con_jacobian!), typeof(Main.var"##226".con_h!), Float64}, initial_x::Vector{Float64}, method::IPNewton{typeof(Optim.backtrack_constrained_grad), Symbol}) @ Optim ~/work/Optim.jl/Optim.jl/src/multivariate/solvers/constrained/ipnewton/interior.jl:229 [4] top-level scope @ ~/work/Optim.jl/Optim.jl/docs/src/examples/generated/ipnewton_basics.ipynb:4 [5] eval @ ./boot.jl:385 [inlined] [6] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String) @ Base ./loading.jl:2076 [7] #44 @ ~/.julia/packages/Literate/kaTtz/src/Literate.jl:861 [inlined] [8] (::IOCapture.var"#4#7"{Core.TypeofBottom, Literate.var"#44#45"{String, Bool, Module, String}, IOContext{Base.PipeEndpoint}, IOContext{Base.PipeEndpoint}, Base.PipeEndpoint, Base.PipeEndpoint})() @ IOCapture ~/.julia/packages/IOCapture/Rzdxd/src/IOCapture.jl:161 [9] with_logstate(f::Function, logstate::Any) @ Base.CoreLogging ./logging.jl:515 [10] with_logger @ ./logging.jl:627 [inlined] [11] capture(f::Literate.var"#44#45"{String, Bool, Module, String}; rethrow::Type, color::Bool, passthrough::Bool, capture_buffer::IOBuffer) @ IOCapture ~/.julia/packages/IOCapture/Rzdxd/src/IOCapture.jl:158 [12] execute_block(sb::Module, block::String; inputfile::String, fake_source::String, softscope::Bool) @ Literate ~/.julia/packages/Literate/kaTtz/src/Literate.jl:859 [13] execute_block @ ~/.julia/packages/Literate/kaTtz/src/Literate.jl:844 [inlined] [14] execute_notebook(nb::Dict{Any, Any}; inputfile::String, fake_source::String, softscope::Bool) @ Literate ~/.julia/packages/Literate/kaTtz/src/Literate.jl:762 [15] (::Literate.var"#38#40"{Dict{String, Any}})() @ Literate ~/.julia/packages/Literate/kaTtz/src/Literate.jl:741 [16] cd(f::Literate.var"#38#40"{Dict{String, Any}}, dir::String) @ Base.Filesystem ./file.jl:112 [17] jupyter_notebook(chunks::Vector{Literate.Chunk}, config::Dict{String, Any}) @ Literate ~/.julia/packages/Literate/kaTtz/src/Literate.jl:740 [18] notebook(inputfile::String, outputdir::String; config::Dict{Any, Any}, kwargs::@Kwargs{execute::Bool}) @ Literate ~/.julia/packages/Literate/kaTtz/src/Literate.jl:677 [19] top-level scope @ ~/work/Optim.jl/Optim.jl/docs/generate.jl:18 [20] include(fname::String) @ Base.MainInclude ./client.jl:489 [21] top-level scope @ ~/work/Optim.jl/Optim.jl/docs/make.jl:10 [22] include(mod::Module, _path::String) @ Base ./Base.jl:495 [23] exec_options(opts::Base.JLOptions) @ Base ./client.jl:318 [24] _start() @ Base ./client.jl:552
* Status: success * Candidate solution Final objective value: 2.966216e-01 * Found with Algorithm: Interior Point Newton * Convergence measures |x - x'| = 0.00e+00 ≤ 0.0e+00 |x - x'|/|x'| = 0.00e+00 ≤ 0.0e+00 |f(x) - f(x')| = 0.00e+00 ≤ 0.0e+00 |f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00 |g(x)| = 7.71e-01 ≰ 1.0e-08 * Work counters Seconds run: 0 (vs limit Inf) Iterations: 34 f(x) calls: 158 ∇f(x) calls: 158
Note that the algorithm warns that the Initial guess is not an
interior point. IPNewton
can often handle this, however, if the
initial guess is such that c(x) = u_c
, then the algorithm currently
fails. We may fix this in the future.
The following example illustrates how to add an additional constraint. In particular, we add a constraint function $$ c(x)_2 = x_2\sin(x_1)-x_1 $$
function con2_c!(c, x)
c[1] = x[1]^2 + x[2]^2 ## First constraint
c[2] = x[2]*sin(x[1])-x[1] ## Second constraint
c
end
function con2_jacobian!(J, x)
# First constraint
J[1,1] = 2*x[1]
J[1,2] = 2*x[2]
# Second constraint
J[2,1] = x[2]*cos(x[1])-1.0
J[2,2] = sin(x[1])
J
end
function con2_h!(h, x, λ)
# First constraint
h[1,1] += λ[1]*2
h[2,2] += λ[1]*2
# Second constraint
h[1,1] += λ[2]*x[2]*-sin(x[1])
h[1,2] += λ[2]*cos(x[1])
# Symmetrize h
h[2,1] = h[1,2]
h
end;
We generate the constraint objects and call IPNewton
with
initial guess $x_0 = (0.25,0.25)$.
x0 = [0.25, 0.25]
lc = [-Inf, 0.0]; uc = [0.5^2, 0.0]
dfc = TwiceDifferentiableConstraints(con2_c!, con2_jacobian!, con2_h!,
lx, ux, lc, uc)
res = optimize(df, dfc, x0, IPNewton())
* Status: success * Candidate solution Final objective value: 1.000000e+00 * Found with Algorithm: Interior Point Newton * Convergence measures |x - x'| = 6.90e-10 ≰ 0.0e+00 |x - x'|/|x'| = 3.55e+08 ≰ 0.0e+00 |f(x) - f(x')| = 1.38e-09 ≰ 0.0e+00 |f(x) - f(x')|/|f(x')| = 1.38e-09 ≰ 0.0e+00 |g(x)| = 2.00e+00 ≰ 1.0e-08 * Work counters Seconds run: 0 (vs limit Inf) Iterations: 29 f(x) calls: 215 ∇f(x) calls: 215
This notebook was generated using Literate.jl.