Worked tutorial on low-level AbstractNLPEvaluator interface?

Specific Domains Optimization (Mathematical)
1

I have a very old lab-internal package that performs a variety of image registration tasks. It was written against MathProgBase, and clearly it’s well past time to rewrite it for MathOptInterface. While the documentation is extensive (thanks!!), I wonder if anyone has a complete tutorial walking through how one solves a problem this way?

I’m struggling even with basics like how to set a regular Julia function as the objective. In case it helps, let me direct you to our old MathProgBase code. As a word of warning, it was never polished for external use, so it’s not exactly easy going. The source file actually has optimizers for several different problems. To pull out the easiest one, there’s an optimizer for rigid-body alignment (rotations + translations). Focusing on just the key points, the two images to be registered are packed into structures as fixed and moving and we’ve written Julia functions to evaluate the objective and to use ForwardDiff to compute derivatives. Finally, here’s the problem setup.

If anyone has an illustration at this level in their own code but using MathOptInterface instead of MathProgBase, I’d be very grateful!

2

If you want to solve an optimization problem with the solvers accessible through MathOptInterface, a good start is trying JuMP first.

You can pass matrices if that is the reason to use MathProgBase: Transitioning from MathProgBase · JuMP

1 Like
3

Without a full example of the code, something like this would work. JuMP can automatically build the ForwardDiff gradient for you. It’s pretty readable what is going on, and you can easily add other constraints.

import JuMP, Ipopt

function solve(d::RigidValue, p0::Vector{Float64}, ub::Vector{Float64})
    @assert length(p0) == length(ub)
    N = length(p0)
    model = JuMP.Model(Ipopt.Optimizer)
    JuMP.@variable(model, -ub[i] <= x[i in 1:N] <= ub[i], start = p0[i])
    JuMP.set_attribute(model, "hessian_approximation", "limited-memory")
    JuMP.set_attribute(model, "sb", "yes")
    JuMP.@operator(model, op_objective, N, (x...) -> d(collect(x)))
    JuMP.@objective(model, Min, op_objective(x...))
    JuMP.optimize!(model)
    if JuMP.termination_status(model) != JuMP.LOCALLY_SOLVED
        @warn("Solution was not optimal")
    end
    p = JuMP.value.(x)
    fval = JuMP.objective_value(model)
    return p, fval
end

but if you don’t have any constraints, then you could even just use Ipopt directly. I didn’t test, but something like this:

function solve(d::RigidValue, p0::Vector{Float64}, ub::Vector{Float64})
    @assert length(p0) == length(ub)
    N = length(p0)
    function eval_grad_f(x::Vector{Float64}, grad_f::Vector{Float64})
        grad_f .= ForwardDiff.gradient(d, x)
        return
    end    
    prob = Ipopt.CreateIpoptProblem(
        N,      # number of variables
        -ub,    # variable lower bound
        ub,     # variable upper bound
        m,      # number of constraints
        Cdouble[],    # constraint lower
        Cdouble[],    # constraint upper
        0,            # non-zeros in Jacobian
        0,            # non-zeros in Hessian
        x::Vector{Cdouble} -> d(x),      # eval objective
        (x, g) -> nothing,               # eval constraints,
        eval_grad_f,                     # eval objective gradient
        (args...) -> nothing,            # eval Jacobian
        nothing,                         # Disable explicit Hessian
    )
    prob.x = p0
    status = Ipopt.IpoptSolve(prob)
    if status != 0
        @warn("Something went wrong")
    end
    return prob.obj_val, prob.x
end
3 Likes
4

Thanks both! I missed the significance of the @operator interface, I think I got convinced that the only option was to use expressions. Once I get it working I will submit a doc PR that would have clarified this issue (and any others I encounter) for me.

EDIT: I think I’ll just use Ipopt directly, no need to answer the question below. Thanks @odow!

With the JuMP-based approach, is there any way I can block it from calling ForwardDiff and instead use my own gradient-evaluation code? I have another “suite” in the same package, focused on GridDeformation, where I compute gradients in a special manner to handle discontinuities in interpolation. I have defined, e.g., MOI.eval_objective_gradient(d::DeformOpt, grad_f, x), and I’ve also defined MOI.features_available(::DeformOpt) = [:Grad, :Jac].

Or perhaps it’s better to just use Ipopt directly? I do have constraints in some problems, but presumably I can elaborate on the example from @odow. What choice would the experts make?

5

With the JuMP-based approach, is there any way I can block it from calling ForwardDiff and instead use my own gradient-evaluation code?

You can provide gradients (and optionally, Hessians): Nonlinear Modeling · JuMP

Or perhaps it’s better to just use Ipopt directly?
What choice would the experts make?

I would encourage almost everybody to use JuMP and try to be clever with hard-coded derivatives etc. The risk of making a mistake with the Jacobians or sparsity is very high. But you are not a novice Julia user…so I don’t have a recommendation.

Really though, if you have Julia functions for the objective, and possibly constraints, and you’re happy constructing your own derivatives, then JuMP is a big heavy package for limited utility. See:

You might be better using something like NLopt.