Description: An example of implementing a derivative-based nonlinear solver in Julia and hooking it into MathProgBase. MathProgBase connects solvers to modeling interfaces like JuMP and AMPL, which provide automatic computation of exact first and second-order derivatives.
Author: Miles Lubin
License:
This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
In this notebook, we demonstrate how to implement Newton's method, by querying derivatives through the MathProgBase nonlinear interface. We then demonstrate using this new solver from both JuMP and AMPL.
using MathProgBase
First we define the NewtonSolver
object, which holds solver options (here there are none). Solver objects are used to instantiate a NewtonData
object, which is the object which then solves the instance of the optimization problem. The NewtonData
type is a subtype of the AbstractNonlinearModel
type and stores all the instance data we need to run the algorithm.
type NewtonSolver <: MathProgBase.AbstractMathProgSolver
end
type NewtonData <: MathProgBase.AbstractNonlinearModel
numVar::Int
d # NLP evaluator
x::Vector{Float64}
hess_I::Vector{Int} # 1st component of Hessian sparsity pattern
hess_J::Vector{Int} # 2nd component of Hessian sparsity pattern
status::Symbol
end
MathProgBase.NonlinearModel(solver::NewtonSolver) = NewtonData(0,nothing,Float64[],Int[],Int[],:Uninitialized);
This method is called when we're about to solve the optimization problem. It provides basic problem dimensions and lower and upper bound vectors, as well as the AbstractNLPEvaluator
object, which is an oracle that can be used by the solver to query function and derivative evaluations.
function MathProgBase.loadproblem!(m::NewtonData, numVar, numConstr, l, u, lb, ub, sense, d::MathProgBase.AbstractNLPEvaluator)
@assert numConstr == 0 # we don't handle constraints
@assert all(l .== -Inf) && all(u .== Inf) # or variable bounds
@assert sense == :Min # or maximization
MathProgBase.initialize(d, [:Grad, :Hess]) # request gradient and hessian evaluations
m.d = d
m.numVar = numVar
m.x = zeros(numVar)
m.hess_I, m.hess_J = MathProgBase.hesslag_structure(d)
end;
And here's the actual implementation of Newton's method. We assume:
We also fix the step size to 1 and do not perform a line search.
function MathProgBase.optimize!(m::NewtonData)
iteration = 0
∇f = Array(Float64,m.numVar)
hess_val = Array(Float64, length(m.hess_I)) # array to store nonzeros in Hessian
MathProgBase.eval_grad_f(m.d, ∇f, m.x) # writes gradient to ∇f vector
∇f_norm = norm(∇f)
while ∇f_norm > 1e-5
println("Iteration $iteration. Gradient norm $∇f_norm, x = $(m.x)");
# Evaluate the Hessian of the Lagrangian.
# We don't have any constraints, so this is just the Hessian.
MathProgBase.eval_hesslag(m.d, hess_val, m.x, 1.0, Float64[])
# The derivative evaluator provides the Hessian matrix
# in lower-triangular sparse (i,j,value) triplet format,
# so now we convert this into a Julia symmetric sparse matrix object.
H = Symmetric(sparse(m.hess_I,m.hess_J,hess_val),:L)
m.x = m.x - H\∇f # newton step
MathProgBase.eval_grad_f(m.d, ∇f, m.x)
∇f_norm = norm(∇f)
iteration += 1
end
println("Iteration $iteration. Gradient norm $∇f_norm, x = $(m.x), DONE");
m.status = :Optimal
end;
We implement a few methods to provide starting solutions and query solution status.
MathProgBase.setwarmstart!(m::NewtonData,x) = (m.x = x)
MathProgBase.status(m::NewtonData) = m.status
MathProgBase.getsolution(m::NewtonData) = m.x
MathProgBase.getconstrduals(m::NewtonData) = [] # lagrange multipliers on constraints (none)
MathProgBase.getreducedcosts(m::NewtonData) = zeros(m.numVar) # lagrange multipliers on variable bounds
MathProgBase.getobjval(m::NewtonData) = MathProgBase.eval_f(m.d, m.x);
using JuMP
jumpmodel = Model(solver=NewtonSolver())
@variable(jumpmodel, x)
@variable(jumpmodel, y)
@NLobjective(jumpmodel, Min, (x-5)^2 + (y-3)^2)
status = solve(jumpmodel);
Iteration 0. Gradient norm 11.661903789690601, x = [0.0,0.0] Iteration 1. Gradient norm 1.9860273225978185e-15, x = [4.999999999999999,2.9999999999999996], DONE
Note that we've converged to the optimal solution in one step. This is expected when applying Newton's method to a strictly convex quadratic objective function.
AMPL is a commercial algebraic modeling language. You can download a limited demo of AMPL here and extract it into the ampl-demo
directory under your home directory to follow along.
mod = """
var x;
var y;
minimize OBJ:
(x-5)^2 + (y-3)^2;
write gquad;
"""
fd = open("quad.mod","w")
write(fd, mod)
close(fd)
The proprietary ampl
executable translates the .mod
file into a .nl
file which we can directly read from Julia.
;~/ampl-demo/ampl quad.mod
;cat quad.nl | head -n 3
The AmplNLReader and AMPLMathProgInterface packages are not yet officially registered as Julia packages. They're not yet recommended for general use.
But anyway, here's a demo of what they can do:
using AmplNLReader, AMPLMathProgInterface
Here we read in the .nl
file and call the AMPLMathProgInterface.loadamplproblem!
method to load the .nl
file into our solver which we just wrote. Then we just call optimize!
and we're done.
nlp = AmplNLReader.AmplModel("quad.nl")
m = MathProgBase.model(NewtonSolver())
AMPLMathProgInterface.loadamplproblem!(m, nlp)
MathProgBase.optimize!(m);
Again we converge to the optimal solution in one step, but this time we used AMPL to compute the derivatives of the model instead of JuMP.