# Worked tutorial on low-level AbstractNLPEvaluator interface?

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!

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

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)
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,
(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

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?

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.