User-defined Hessians

This tutorial was generated using Literate.jl. Download the source as a .jl file.

In this tutorial, we explain how to write a user-defined operator (see User-defined operators) with a Hessian matrix explicitly provided by the user.

For a more advanced example, see Nested optimization problems.

Required packages

This tutorial uses the following packages:

using JuMP
import Ipopt

Rosenbrock example

As a simple example, we consider the Rosenbrock function:

rosenbrock(x...) = (1 - x[1])^2 + 100 * (x[2] - x[1]^2)^2

rosenbrock (generic function with 1 method)

which has the gradient vector:

function ∇rosenbrock(g::AbstractVector, x...)
    g[1] = 400 * x[1]^3 - 400 * x[1] * x[2] + 2 * x[1] - 2
    g[2] = 200 * (x[2] - x[1]^2)
    return
end

∇rosenbrock (generic function with 1 method)

and the Hessian matrix:

function ∇²rosenbrock(H::AbstractMatrix, x...)
    H[1, 1] = 1200 * x[1]^2 - 400 * x[2] + 2
    # H[1, 2] = -400 * x[1] <-- not needed because Hessian is symmetric
    H[2, 1] = -400 * x[1]
    H[2, 2] = 200.0
    return
end

∇²rosenbrock (generic function with 1 method)

You may assume the Hessian matrix H is initialized with zeros, and because it is symmetric you need only to fill in the non-zero of the lower-triangular terms.

The matrix type passed in as H depends on the automatic differentiation system, so make sure the first argument to the Hessian function supports an AbstractMatrix (it may be something other than Matrix{Float64}). However, you may assume only that H supports size(H) and setindex!.

Now that we have the function, its gradient, and its Hessian, we can construct a JuMP model, add the operator, and use it in a macro:

model = Model(Ipopt.Optimizer)
@variable(model, x[1:2])
@operator(model, op_rosenbrock, 2, rosenbrock, ∇rosenbrock, ∇²rosenbrock)
@objective(model, Min, op_rosenbrock(x[1], x[2]))
optimize!(model)
assert_is_solved_and_feasible(model)
solution_summary(model; verbose = true)

solution_summary(; result = 1, verbose = true)
├ solver_name          : Ipopt
├ Termination
│ ├ termination_status : LOCALLY_SOLVED
│ ├ result_count       : 1
│ └ raw_status         : Solve_Succeeded
├ Solution (result = 1)
│ ├ primal_status        : FEASIBLE_POINT
│ ├ dual_status          : FEASIBLE_POINT
│ ├ objective_value      : 1.21190e-27
│ ├ dual_objective_value : 0.00000e+00
│ └ value
│   ├ x[1] : 1.00000e+00
│   └ x[2] : 1.00000e+00
└ Work counters
  ├ solve_time (sec)   : 4.08580e-02
  └ barrier_iterations : 14